#' ---
#' title: "レプリケーション：福祉国家に対する態度決定要因としての社会保障の普遍性と消費税の逆進性"
#' author: "加藤言人"
#' date: "2021年7月28日"
#' number_sections: true
#' toc: true
#' documentclass: bxjsarticle
#' classoption: a4paper, xelatex, ja=standard
#' mainfont: IPAGothic
#' monofont: IPAGothic
#' ---

#+ warning=FALSE, echo=FALSE
# このドキュメントはUTF-8で書かれています。
# This document is encoded in UTF-8

#+ warning=FALSE, echo=FALSE
## Set Working Directory to the current directory 
## (If using RStudio, can be set automatically)
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

## Clear Workspace
rm(list=ls())

## Used Packages
library(estvis); library(texreg)
library(readr);library(ggplot2)
require(grid);library(gridExtra);library(haven)
require(sandwich); require(lmtest); require(xtable)

library(knitr)
# 出力フォーマットが TeX（PDF含む）の場合のみ対処する
if (knitr::opts_knit$get("rmarkdown.pandoc.to") %in% c("beamer", "latex")) {
  
  # conversion failure on '...' in 'mbcsToSbcs' の Warning 発生の workaround
  options(device = function(file, width = 7, height = 7, ...) {
    cairo_pdf(tempfile(), width = width, height = height, ...)
  })
  
  ## 1. cairo_pdf を使う方法
  # * family には OS にインストールされているフォント名を指定する。
  knitr::opts_chunk$set(dev="cairo_pdf", dev.args=list(family="IPAGothic"))
}

### Font Settings
# windowsFonts(YuGothic = windowsFont("Yu Gothic"))
# theme_set(theme_gray(base_size = 10, base_family = "YuGothic"))

#+ warning=FALSE, message=FALSE
## Read Data
d <- readRDS("main_data_tax_v3.rds")

## Drop Satisficers
d <- d[which(d$satisficer==0),] 

## Keep Respondents without missing values for relevant variables
dtmp <- d[complete.cases(d[,c("tax1_opi","tax2_opi","inc",
                              "knall","fem","age","lvlen","ownh",
                              "edu3","wk","mar","cld")]),]
# nrow(dtmp)

#'
#' # イデオロギー指標間の相関（表A2）
#' 

ctab <- cor(dtmp[,c("ide_self","ide_psup","ide_iss_1","ide_iss_2")])
ctab[upper.tri(ctab)] <- NA
colnames(ctab) <- rownames(ctab) <- c("自己申告","政党支持","外交安全保障","権利機会平等")

## 表A2
round(ctab, 3)
# print(xtable(ctab, digits=3,caption="イデオロギー指標間の相関"), 
#       caption.placement="top")


#'
#' # 記述統計
#' 
#' ## 従属変数
#'
#' ### 理想消費税率の分布（図A1）
#' 

tab <- data.frame(tax_opi=c(dtmp$tax1_opi,rep(NA, nrow(dtmp)),
                            dtmp$tax2_opi,rep(NA, nrow(dtmp))),
                  tax_opi_sq=c(rep(NA, nrow(dtmp)),sqrt(dtmp$tax1_opi),
                               rep(NA, nrow(dtmp)),sqrt(dtmp$tax2_opi)),
                  taxcat=rep(c("食料・日用品","他全ての商品"),
                             each=nrow(dtmp)*2),
                  scale=rep(c("スケール未変換","平方根スケール"),
                            each=nrow(dtmp)))
tab$taxcat <- factor(tab$taxcat, levels=unique(tab$taxcat))

tab2 <- data.frame(tax_opi_m = c(mean(dtmp$tax1_opi, na.rm=TRUE),
                                 mean(sqrt(dtmp$tax1_opi), na.rm=TRUE),
                                 mean(dtmp$tax2_opi, na.rm=TRUE),
                                 mean(sqrt(dtmp$tax2_opi), na.rm=TRUE)),
                   taxcat=rep(c("食料・日用品","他全ての商品"),
                              each=2),
                   scale=rep(c("スケール未変換","平方根スケール"),
                             2),
                   y = c(0.4,0.02,0.4,0.02))
tab2$taxcat <- factor(tab2$taxcat, levels=unique(tab$taxcat))
tab2$tax_opi_lab <- paste0("平均: ", sprintf("%2.2f",tab2$tax_opi_m))

## 図A1
p <- ggplot(tab) + 
  geom_histogram(aes(x=tax_opi, y=(..count../sum(..count..))*2), bins=11, alpha=0.5) +
  geom_density(aes(x=tax_opi, y=..density..*4), bw=2) +
  geom_histogram(aes(x=tax_opi_sq, y=(..count../sum(..count..))*2), bins=11, alpha=0.5) +
  geom_density(aes(x=tax_opi_sq), bw=0.5) +
  geom_vline(data=tab2, aes(xintercept=tax_opi_m), color="red", linetype=2) + 
  geom_text(data=tab2, aes(label=tax_opi_lab, x=tax_opi_m, y=y), hjust=-0.05) + 
  facet_grid(taxcat~scale, scales = "free_x") + 
  scale_y_continuous(breaks = c(0,0.1,0.2,0.3,0.4), 
                     labels = c(0,0.05,0.1,0.15,0.2)) + 
  # scale_x_continuous(breaks=sqrt(c(0,2.5,5,10,25,50,100)),
  #                    labels=c(0,1,5,10,25,50,100)) + 
  ylab("回答割合") + xlab("理想消費税率\n（軸のスケール：左＝そのままの値、右＝平方根）") + 
  #ggtitle("理想消費税率の回答分布") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust=0.5, face="bold"),
        axis.text.x = element_text(size=12, face="bold"))

#+ warning=FALSE, fig.width=7, fig.height=5
p

#+ warning=FALSE, eval=FALSE
# ggsave("dvdist.png", p, width=7, height=5)


#'
#' ## 独立変数：収入とイデオロギー
#' 
#' ### 収入
#'

tab <- table(dtmp$inc)/sum(table(dtmp$inc))
tab <- data.frame(prop = as.numeric(tab),
                  names = c("～\n200万円\n(1)",
                            "200〜\n400万円\n(2)",
                            "400～\n600万円\n(3)",
                            "600〜\n800万円\n(4)",
                            "800〜\n1000万円\n(5)",
                            "1000〜\n1200万円\n(6)",
                            "1200〜\n1400万円\n(7)",
                            "1400万円\n～\n(8)"))
tab$names <- factor(tab$names, levels=tab$names)

## 図A2（上側）
p0 <- ggplot(tab, aes(x=names,y=prop)) + 
  geom_bar(stat="identity") + 
  ylab(NULL) + xlab(NULL) + 
  ggtitle("世帯収入（度数分布）") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust=0.5, face="bold"),
        axis.text.x = element_text(size=12, face="bold"))
# p0

#' 
#' ### 自己申告イデオロギー
#'
tab <- table(dtmp$ide_self)/sum(table(dtmp$ide_self))
tab <- data.frame(prop = as.numeric(tab),
                  names = c("左派/\nリベラル\n(-3)","-2","-1",
                            "中立\n(0)","1","2","右派/\n保守\n(3)"))
tab$names <- factor(tab$names, levels=tab$names)

## 図A2（下側）
p1 <- ggplot(tab, aes(x=names,y=prop)) + 
  geom_bar(stat="identity") + 
  ylab(NULL) + xlab(NULL) + 
  ggtitle("自己申告イデオロギー（度数分布）") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust=0.5, face="bold"),
        axis.text.x = element_text(size=12, face="bold"))
# p1

#'
#' ### 世帯収入と自己申告イデオロギーの分布（図A2）
#'

#+ warning=FALSE, fig.width=8, fig.height=6
ggplot() + theme_void()
p <- arrangeGrob(p0,p1, nrow=2, left="回答割合")
grid.draw(p)

#+ warning=FALSE, eval=FALSE
# ggsave("incide.png", p, width=8, height=6)

#' 
#' ### 争点態度イデオロギー
#' 

## 図A7（中央）
p2_1 <- ggplot(dtmp[,"ide_iss_1"],
               aes(x=ide_iss_1,y=..count../sum(..count..))) +
  geom_histogram(bins=10, color="white") +
  ylab(NULL) + xlab(NULL) + 
  ggtitle("外交安全保障\nイデオロギー\n（ヒストグラム）") + 
  scale_x_continuous(breaks=c(-3,-2,-1,0,1,2,3),
                     limits=c(-3,3),
                     labels=c("左派\n(-3)\n","-2","-1","0","1","2","右派\n(3)\n")) +
  theme_bw() + 
  theme(plot.title = element_text(hjust=0.5, face="bold"),
        axis.text.x = element_text(size=12, face="bold"))
# p2_1

## 図A7（右）
p2_2 <- ggplot(dtmp[,"ide_iss_2"],
               aes(x=ide_iss_2,y=..count../sum(..count..))) + 
  geom_histogram(bins=10,color="white") +
  ylab(NULL) + xlab(NULL) + 
  ggtitle("権利機会平等\nイデオロギー\n（ヒストグラム）") + 
  scale_x_continuous(breaks=c(-3,-2,-1,0,1,2,3),
                     limits=c(-3,3),
                     labels=c("左派\n(-3)\n","-2","-1","0","1","2","右派\n(3)\n")) +
  theme_bw() + 
  theme(plot.title = element_text(hjust=0.5, face="bold"),
        axis.text.x = element_text(size=12, face="bold"))
# p2_2

#'
#' ### 政党支持イデオロギー
#' 

## 図A7（左）
tab <- table(dtmp$ide_psup)/sum(table(dtmp$ide_psup))
tab <- data.frame(prop = as.numeric(tab),
                  names = c("左派\n政党支持\n(-1)","無党派\nその他\n(0)",
                            "右派\n政党支持\n(1)"))
tab$names <- factor(tab$names, levels=tab$names)

p3 <- ggplot(tab, aes(x=names,y=prop)) + 
  geom_bar(stat="identity") + 
  ylab(NULL) + xlab(NULL) + 
  ggtitle("政党支持\nイデオロギー\n（度数分布）") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust=0.5, face="bold"),
        axis.text.x = element_text(size=12, face="bold"))
# p3

#'
#' ### 政党支持・争点態度イデオロギー変数の分布（図A7）
#' 

#+ warning=FALSE, fig.width=8.5, fig.height=4
ggplot() + theme_void()
p <- arrangeGrob(p3,p2_1,p2_2, nrow=1, left="回答割合")
grid.draw(p)

#+ warning=FALSE, eval=FALSE
# ggsave("idedist32.png", p, width=8.5, height=4)

#'
#' # 実験群比較
#'
#' ## 準備

# 統制変数
ctl <- formula( ~ .+ knall + fem + age + lvlen + ownh + 
                  as.factor(edu3) + wk + mar + cld)


# 変数名
vn <- c("（定数項）",
        "1.逆進性",
        "2.社会保障普遍性",
        "3.社会保障選別性", 
        "4.逆進性＆社会保障普遍性",
        "5.逆進性＆社会保障選別性",
        "政治知識","性別（女性）",
        "年齢","居住年数","持ち家",
        "教育：短大／高専／専門学校",
        "教育：大卒以上",
        "就労","婚姻","子ども")

vnx <- c(vn[1:6],"イデオロギー",vn[7:16],
         "イデオロギー×1.逆進",
         "イデオロギー×2.普遍",
         "イデオロギー×3.選別",
         "イデオロギー×4.逆進＆普遍",
         "イデオロギー×5.逆進＆選別",
         "イデオロギー",
         "イデオロギー×1.逆進",
         "イデオロギー×2.普遍",
         "イデオロギー×3.選別",
         "イデオロギー×4.逆進＆普遍",
         "イデオロギー×5.逆進＆選別",
         "イデオロギー",
         "イデオロギー×1.逆進",
         "イデオロギー×2.普遍",
         "イデオロギー×3.選別",
         "イデオロギー×4.逆進＆普遍",
         "イデオロギー×5.逆進＆選別")
vnx2 <- c(vn,vnx[c(7,18:34)])

vnxinc <- c(vn[1:6],"世帯収入",vn[7:16],
            "世帯収入×1.逆進",
            "世帯収入×2.普遍", 
            "世帯収入×3.選別",
            "世帯収入×4.逆進＆普遍",
            "世帯収入×5.逆進＆選別")

#'
#' ## バランスチェック（図A3）
#'

## バランスチェック関数
checkbal <- function(dtlist, dtnames, robust=TRUE){
  
  restab <- list()
  datasub <- list()
  for(i in 1:length(dtlist)){
    
    dt <- dtlist[[i]]
    
    dmod <- data.frame(
      knall = dt$knall,
      fem0 = ifelse(dt$fem==0,1,0),
      fem1 = ifelse(dt$fem==1,1,0),
      ageby10 = dt$age/10,
      lvlen = dt$lvlen,
      ownh0 = ifelse(dt$ownh==0,1,0),
      ownh1 = ifelse(dt$ownh==1,1,0),
      edu30 = ifelse(dt$edu3==0,1,0),
      edu31 = ifelse(dt$edu3==1,1,0),
      edu32 = ifelse(dt$edu3==2,1,0),
      wk0 = ifelse(dt$wk==0,1,0),
      wk1 = ifelse(dt$wk==1,1,0),
      mar0 = ifelse(dt$mar==0,1,0),
      mar1 = ifelse(dt$mar==1,1,0),
      cld0 = ifelse(dt$cld==0,1,0),
      cld1 = ifelse(dt$cld==1,1,0),
      inc = dt$inc,
      ide_self = dt$ide_self,
      ide_psup = dt$ide_psup,
      ide_iss_1 = dt$ide_iss_1,
      ide_iss_2 = dt$ide_iss_2
    )
    
    if (robust==TRUE){
      res1 <- apply(dmod, 2, function(x) coeftest(lm(x~dt$treated),
                                                  vcov.=vcovHC(lm(x~dt$treated),type="HC2")))
    } else {
      res1 <- apply(dmod, 2, function(x) coeftest(lm(x~dt$treated)))
    }
    
    res2 <- as.data.frame(t(apply(res1,2,function(x) c(x[2],x[4],x[8]))))
    colnames(res2) <- c("coef","se","pval")
    res2$lCI <- res2$coef - 1.96*res2$se
    res2$uCI <- res2$coef + 1.96*res2$se
    
    res3 <- data.frame(Variable=
                         rep(
                           c("政治知識",
                             "性別（男性）",
                             "性別（女性）",
                             "年齢（10歳ごと）",
                             "居住年数",
                             "持ち家なし",
                             "持ち家あり",
                             "教育：高卒以下",
                             "教育：短大／高専／専門学校",
                             "教育：大卒以上",
                             "無職",
                             "有職",
                             "結婚していない",
                             "結婚している",
                             "子どもなし",
                             "子どもあり",
                             "世帯収入",
                             "自己申告イデオロギー",
                             "政党支持イデオロギー",
                             "外交安全保障イデオロギー",
                             "権利機会平等イデオロギー"
                           ),2),
                       stat = c(res2$coef,res2$pval),
                       lCI = c(res2$lCI,rep(NA,nrow(res2))),
                       uCI = c(res2$uCI,rep(NA,nrow(res2))),
                       val = rep(c("実験群ー統制群","p値"),each=21)
    )
    res3$Variable <- factor(res3$Variable,levels=rev(unique(res3$Variable)))
    res3$val <- factor(res3$val,levels=unique(res3$val))
    restab[[i]] <- res3
    
  }
  
  if(length(dtlist)==1){
    restab <- restab[[1]]
  } else {
    restabtemp <- restab[[1]]
    restabtemp$group <- dtnames[1]
    for (i in 2:length(dtlist)){
      restabtemp2 <- restab[[i]]
      restabtemp2$group <- dtnames[i]
      restabtemp <- rbind(restabtemp,restabtemp2)
    }
    restab <- restabtemp
  }
  restab$group <- factor(restab$group, unique(restab$group))
  
  data2 <- data.frame(val=c("実験群ー統制群","p値"),
                      vloc1=c(NA,0),
                      vloc2=c(0,0.1))
  data2$val <- factor(data2$val,levels=unique(data2$val))
  
  p <- ggplot(restab, aes(x=Variable,y=stat)) + 
    geom_errorbar(aes(ymin=lCI,ymax=uCI, color=group),width=0.5, 
                  position=position_dodge(width=-0.7)) +
    geom_point(aes(shape=group, color=group), 
               position=position_dodge(width=-0.7)) + 
    geom_hline(data=data2, aes(yintercept=vloc1), linetype=1) +
    geom_hline(data=data2, aes(yintercept=vloc2), linetype=2) +
    scale_y_continuous(breaks=c(-0.1,0,0.1,0.3,0.6,0.9)) +
    scale_shape_discrete(name="実験群") + 
    scale_color_brewer(name="実験群", palette=2, type = "qual") +
    facet_grid(.~val, scales = "free_x") + coord_flip() + 
    xlab(NULL) + ylab(NULL) + 
    labs(caption = "注：各変数を従属変数、実験群を従属変数としてOLS回帰し、各実験群と統制群の差を推定。横線は95%信頼区間（ロバスト標準誤差を使用）。") +
    theme_bw() + 
    theme(legend.position = "bottom") 
  
  return(p)
}

dlist <- lapply(c(1,2,3,4,5), function(k) {
  out <- subset(d, g_ctax_N%in%c(0,k))
  out$treated = ifelse(out$g_ctax_N>0,1,0)
  return(out)
})
dnames <- c("1.逆進性","2.社会保障普遍性","3.社会保障選別性",
            "4.逆進性&社会保障普遍性","5.逆進性&社会保障選別性")

# 図A3
pbalance <- checkbal(dlist,dnames)

#+ warning=FALSE, fig.width=10, fig.height=10
pbalance

#+ warning=FALSE, eval=FALSE, echo=FALSE
# ggsave("balcheck.png", pbalance, width=10, height=10)

#'
#' ## 実験刺激の直接効果
#'

# 推定
m_ctax1 <- lm(update(tax1_opi ~ as.factor(g_ctax_N),ctl), data=dtmp)
m0_1 <- m_ctax1
m_ctax2 <- lm(update(tax2_opi ~ as.factor(g_ctax_N),ctl), data=dtmp)
m0_2 <- m_ctax2

coeftest(m_ctax1, vcov.=vcovHC(m_ctax1,"HC2"))
coeftest(m_ctax2, vcov.=vcovHC(m_ctax2,"HC2"))

# 普遍性刺激の効果（４種類）
# 統制群vs実験群2＆実験群3vs実験群2＆実験群1vs実験群4＆実験群5vs実験群4
m0_h10_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5)),ctl), data=dtmp)
m0_h10_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5)),ctl), data=dtmp)
m0_h11_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5)),ctl), data=dtmp)
m0_h11_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5)),ctl), data=dtmp)
m0_h12_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5)),ctl), data=dtmp)
m0_h12_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5)),ctl), data=dtmp)
m0_h13_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3)),ctl), data=dtmp)
m0_h13_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3)),ctl), data=dtmp)

m0r_h10_1 <- coeftest(m0_h10_1, vcov.=vcovHC(m0_h10_1,"HC2"))
m0r_h10_2 <- coeftest(m0_h10_2, vcov.=vcovHC(m0_h10_2,"HC2"))
m0r_h11_1 <- coeftest(m0_h11_1, vcov.=vcovHC(m0_h11_1,"HC2"))
m0r_h11_2 <- coeftest(m0_h11_2, vcov.=vcovHC(m0_h11_2,"HC2"))
m0r_h12_1 <- coeftest(m0_h12_1, vcov.=vcovHC(m0_h12_1,"HC2"))
m0r_h12_2 <- coeftest(m0_h12_2, vcov.=vcovHC(m0_h12_2,"HC2"))
m0r_h13_1 <- coeftest(m0_h13_1, vcov.=vcovHC(m0_h13_1,"HC2"))
m0r_h13_2 <- coeftest(m0_h13_2, vcov.=vcovHC(m0_h13_2,"HC2"))

# 逆進性刺激の効果（3種類）
# 統制群vs実験群1＆実験群2vs実験群4 ＆ 実験群3vs実験群5
m0_h20_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5)),ctl), data=dtmp)
m0_h20_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5)),ctl), data=dtmp)
m0_h21_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5)),ctl), data=dtmp)
m0_h21_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5)),ctl), data=dtmp)
m0_h22_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4)),ctl), data=dtmp)
m0_h22_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4)),ctl), data=dtmp)

m0r_h20_1 <- coeftest(m0_h20_1, vcov.=vcovHC(m0_h20_1,"HC2"))
m0r_h20_2 <- coeftest(m0_h20_2, vcov.=vcovHC(m0_h20_2,"HC2"))
m0r_h21_1 <- coeftest(m0_h21_1, vcov.=vcovHC(m0_h21_1,"HC2"))
m0r_h21_2 <- coeftest(m0_h21_2, vcov.=vcovHC(m0_h21_2,"HC2"))
m0r_h22_1 <- coeftest(m0_h22_1, vcov.=vcovHC(m0_h22_1,"HC2"))
m0r_h22_2 <- coeftest(m0_h22_2, vcov.=vcovHC(m0_h22_2,"HC2"))

#' 
#' ### 仮説に関係する実験群比較に関連する直接効果（図A4）
#'

## 作図用データ
htest <- data.frame(dv = rep(c("生活必需品","その他すべて"), each=7),
                    h = rep(c("社会保障普遍性","社会保障普遍性",
                              "社会保障普遍性","社会保障普遍性",
                              "消費税逆進性","消費税逆進性",
                              "消費税逆進性"),2),
                    cp = rep(c("2.普遍 - 0.統制",
                               "2.普遍 - 3.選別",
                               "4.逆進+普遍 - 1.逆進",
                               "4.逆進+普遍 - 5.逆進+選別",
                               "1.逆進 - 0.統制",
                               "4.普遍+逆進 - 2.普遍",
                               "5.選別+逆進 - 3.選別"),2),
                    rbind(m0r_h10_1[2,],m0r_h11_1[2,],m0r_h12_1[2,],m0r_h13_1[2,],
                          m0r_h20_1[2,],m0r_h21_1[2,],m0r_h22_1[2,],
                          m0r_h10_2[2,],m0r_h11_2[2,],m0r_h12_2[2,],m0r_h13_2[2,],
                          m0r_h20_2[2,],m0r_h21_2[2,],m0r_h22_2[2,]))
htest$dv <- factor(htest$dv, levels=unique(htest$dv))
htest$cp <- factor(htest$cp, levels=rev(unique(htest$cp)))
htest$h <- factor(htest$h, levels=unique(htest$h))
htest$lo95 <- htest$Estimate - qnorm(0.975)*htest$Std..Error
htest$up95 <- htest$Estimate + qnorm(0.975)*htest$Std..Error
htest$lo90 <- htest$Estimate - qnorm(0.95)*htest$Std..Error
htest$up90 <- htest$Estimate + qnorm(0.95)*htest$Std..Error

## プロット
p <- ggplot(htest, aes(x=cp, y=Estimate)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lo95,ymax=up95), width=0.25) + 
  geom_errorbar(aes(ymin=lo90,ymax=up90), width=0, size=1.5) + 
  geom_point(size=3) +
  facet_grid(h~dv, scale="free_y", switch="y") + 
  coord_flip() + 
  # scale_color_brewer(name="対象商品", type="qual", palette=2) + 
  labs(x=NULL, y="実験刺激効果（従属変数は理想消費税率）", 
       caption="分析の詳細は回帰表を参照。統制変数有。太線は90%信頼区間、細線は95%信頼区間を示している。", 
       subtitle = "対象商品") + 
  theme_bw() + 
  theme(plot.subtitle = element_text(hjust=0.5),
        strip.background.y = element_blank(),
        strip.text.y = element_blank(),
        strip.placement = "outside")

#+ warning=FALSE, fig.width=8, fig.height=5, dev="png", dpi=300, echo=FALSE
p

#+ warning=FALSE, eval=FALSE
# ggsave("htest_m0_originalscale.png", p, width=8, height=5)

#'
#' ## 世帯収入による条件付け
#'

# 推定
mx_ctax1 <- lm(update(tax1_opi ~ as.factor(g_ctax_N)*inc,ctl), data=dtmp)
coeftest(mx_ctax1, vcovHC(mx_ctax1, "HC2"))
mA_1 <- mx_ctax1
mx_ctax2 <- lm(update(tax2_opi ~ as.factor(g_ctax_N)*inc,ctl), data=dtmp)
coeftest(mx_ctax2, vcovHC(mx_ctax2, "HC2"))
mA_2 <- mx_ctax2

# 仮説１の検証（４種類）
# 統制群vs実験群2＆実験群3vs実験群2＆実験群1vs実験群4＆実験群5vs実験群4
mA1_h10_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(inc-2),ctl), data=dtmp)
mA1_h10_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(inc-2),ctl), data=dtmp)
mA1_h11_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(inc-2),ctl), data=dtmp)
mA1_h11_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(inc-2),ctl), data=dtmp)
mA1_h12_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(inc-2),ctl), data=dtmp)
mA1_h12_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(inc-2),ctl), data=dtmp)
mA1_h13_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(inc-2),ctl), data=dtmp)
mA1_h13_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(inc-2),ctl), data=dtmp)

mA1r_h10_1 <- coeftest(mA1_h10_1, vcov.=vcovHC(mA1_h10_1,"HC2"))
mA1r_h10_2 <- coeftest(mA1_h10_2, vcov.=vcovHC(mA1_h10_2,"HC2"))
mA1r_h11_1 <- coeftest(mA1_h11_1, vcov.=vcovHC(mA1_h11_1,"HC2"))
mA1r_h11_2 <- coeftest(mA1_h11_2, vcov.=vcovHC(mA1_h11_2,"HC2"))
mA1r_h12_1 <- coeftest(mA1_h12_1, vcov.=vcovHC(mA1_h12_1,"HC2"))
mA1r_h12_2 <- coeftest(mA1_h12_2, vcov.=vcovHC(mA1_h12_2,"HC2"))
mA1r_h13_1 <- coeftest(mA1_h13_1, vcov.=vcovHC(mA1_h13_1,"HC2"))
mA1r_h13_2 <- coeftest(mA1_h13_2, vcov.=vcovHC(mA1_h13_2,"HC2"))

mA2_h10_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(inc-5),ctl), data=dtmp)
mA2_h10_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(inc-5),ctl), data=dtmp)
mA2_h11_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(inc-5),ctl), data=dtmp)
mA2_h11_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(inc-5),ctl), data=dtmp)
mA2_h12_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(inc-5),ctl), data=dtmp)
mA2_h12_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(inc-5),ctl), data=dtmp)
mA2_h13_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(inc-5),ctl), data=dtmp)
mA2_h13_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(inc-5),ctl), data=dtmp)

mA2r_h10_1 <- coeftest(mA2_h10_1, vcov.=vcovHC(mA2_h10_1,"HC2"))
mA2r_h10_2 <- coeftest(mA2_h10_2, vcov.=vcovHC(mA2_h10_2,"HC2"))
mA2r_h11_1 <- coeftest(mA2_h11_1, vcov.=vcovHC(mA2_h11_1,"HC2"))
mA2r_h11_2 <- coeftest(mA2_h11_2, vcov.=vcovHC(mA2_h11_2,"HC2"))
mA2r_h12_1 <- coeftest(mA2_h12_1, vcov.=vcovHC(mA2_h12_1,"HC2"))
mA2r_h12_2 <- coeftest(mA2_h12_2, vcov.=vcovHC(mA2_h12_2,"HC2"))
mA2r_h13_1 <- coeftest(mA2_h13_1, vcov.=vcovHC(mA2_h13_1,"HC2"))
mA2r_h13_2 <- coeftest(mA2_h13_2, vcov.=vcovHC(mA2_h13_2,"HC2"))

# 仮説２の検証（３種類）
# 統制群vs実験群1＆実験群2vs実験群4 ＆ 実験群3vs実験群5
mA1_h20_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(inc-2),ctl), data=dtmp)
mA1_h20_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(inc-2),ctl), data=dtmp)
mA1_h21_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(inc-2),ctl), data=dtmp)
mA1_h21_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(inc-2),ctl), data=dtmp)
mA1_h22_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(inc-2),ctl), data=dtmp)
mA1_h22_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(inc-2),ctl), data=dtmp)

mA1r_h20_1 <- coeftest(mA1_h20_1, vcov.=vcovHC(mA1_h20_1,"HC2"))
mA1r_h20_2 <- coeftest(mA1_h20_2, vcov.=vcovHC(mA1_h20_2,"HC2"))
mA1r_h21_1 <- coeftest(mA1_h21_1, vcov.=vcovHC(mA1_h21_1,"HC2"))
mA1r_h21_2 <- coeftest(mA1_h21_2, vcov.=vcovHC(mA1_h21_2,"HC2"))
mA1r_h22_1 <- coeftest(mA1_h22_1, vcov.=vcovHC(mA1_h22_1,"HC2"))
mA1r_h22_2 <- coeftest(mA1_h22_2, vcov.=vcovHC(mA1_h22_2,"HC2"))

mA2_h20_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(inc-5),ctl), data=dtmp)
mA2_h20_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(inc-5),ctl), data=dtmp)
mA2_h21_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(inc-5),ctl), data=dtmp)
mA2_h21_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(inc-5),ctl), data=dtmp)
mA2_h22_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(inc-5),ctl), data=dtmp)
mA2_h22_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(inc-5),ctl), data=dtmp)

mA2r_h20_1 <- coeftest(mA2_h20_1, vcov.=vcovHC(mA2_h20_1,"HC2"))
mA2r_h20_2 <- coeftest(mA2_h20_2, vcov.=vcovHC(mA2_h20_2,"HC2"))
mA2r_h21_1 <- coeftest(mA2_h21_1, vcov.=vcovHC(mA2_h21_1,"HC2"))
mA2r_h21_2 <- coeftest(mA2_h21_2, vcov.=vcovHC(mA2_h21_2,"HC2"))
mA2r_h22_1 <- coeftest(mA2_h22_1, vcov.=vcovHC(mA2_h22_1,"HC2"))
mA2r_h22_2 <- coeftest(mA2_h22_2, vcov.=vcovHC(mA2_h22_2,"HC2"))

#'
#' ### 世帯収入に条件付けされた実験情報刺激の限界効果を用いた仮説検証（図2）
#'

htest <- data.frame(int = rep(c("200〜400万円(10%=2)","800〜1000万円(90%=5)"), each=14), 
                    dv = rep(c("生活必需品","その他すべて"), each=7),
                    h = rep(c("H1A","H1A",
                              "H1A","H1A",
                              "H2A","H2A",
                              "H2A"),4),
                    cp = rep(c("2.普遍 - 0.統制",
                               "2.普遍 - 3.選別",
                               "4.逆進+普遍 - 1.逆進",
                               "4.逆進+普遍 - 5.逆進+選別",
                               "1.逆進 - 0.統制",
                               "4.普遍+逆進 - 2.普遍",
                               "5.選別+逆進 - 3.選別"),4),
                    rbind(mA1r_h10_1[2,],mA1r_h11_1[2,],mA1r_h12_1[2,],mA1r_h13_1[2,],
                          mA1r_h20_1[2,],mA1r_h21_1[2,],mA1r_h22_1[2,],
                          mA1r_h10_2[2,],mA1r_h11_2[2,],mA1r_h12_2[2,],mA1r_h13_2[2,],
                          mA1r_h20_2[2,],mA1r_h21_2[2,],mA1r_h22_2[2,],
                          mA2r_h10_1[2,],mA2r_h11_1[2,],mA2r_h12_1[2,],mA2r_h13_1[2,],
                          mA2r_h20_1[2,],mA2r_h21_1[2,],mA2r_h22_1[2,],
                          mA2r_h10_2[2,],mA2r_h11_2[2,],mA2r_h12_2[2,],mA2r_h13_2[2,],
                          mA2r_h20_2[2,],mA2r_h21_2[2,],mA2r_h22_2[2,]))
htest$dv <- factor(htest$dv, levels=unique(htest$dv))
htest$cp <- factor(htest$cp, levels=rev(unique(htest$cp)))
htest$lo95 <- htest$Estimate - qnorm(0.975)*htest$Std..Error
htest$up95 <- htest$Estimate + qnorm(0.975)*htest$Std..Error
htest$lo90 <- htest$Estimate - qnorm(0.95)*htest$Std..Error
htest$up90 <- htest$Estimate + qnorm(0.95)*htest$Std..Error

p <- ggplot(htest, aes(x=cp, y=Estimate)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lo95,ymax=up95, color=int), width=0.25, position = position_dodge(width=-0.5)) + 
  geom_errorbar(aes(ymin=lo90,ymax=up90, color=int), width=0, size=1.5, position = position_dodge(width=-0.5)) + 
  geom_point(aes(color=int, shape=int), size=3, position = position_dodge(width=-0.5)) +
  facet_grid(h~dv, scales = "free_y", space = "free_y") + 
  coord_flip() + 
  scale_color_brewer(name="世帯収入", type="qual", palette=2) + 
  scale_shape_discrete(name="世帯収入") + 
  labs(x=NULL, y="実験刺激効果（従属変数は理想消費税率）", 
       caption="分析の詳細は回帰表を参照。太線は90%信頼区間、細線は95%信頼区間を示している。",
       subtitle = "対象商品") + 
  theme_bw() + theme(legend.position="bottom",
                     plot.subtitle = element_text(hjust=0.5),
                     strip.placement = "outside")

#+ warning=FALSE, fig.width=8, fig.height=5, dev="png", dpi=300, echo=FALSE
p

#+ warning=FALSE, eval=FALSE
# ggsave("htest_mA_originalscale.png", p, width=8, height=6)

#'
#' ## 自己申告イデオロギーによる条件付け
#'

# 推定
mx_ctax1 <- lm(update(tax1_opi ~ as.factor(g_ctax_N)*ide_self,ctl), data=dtmp)
coeftest(mx_ctax1, vcovHC(mx_ctax1, "HC2"))
mB_1 <- mx_ctax1
mx_ctax2 <- lm(update(tax2_opi ~ as.factor(g_ctax_N)*ide_self,ctl), data=dtmp)
coeftest(mx_ctax2, vcovHC(mx_ctax2, "HC2"))
mB_2 <- mx_ctax2

# 仮説１の検証（４種類）
# 統制群vs実験群2＆実験群3vs実験群2＆実験群1vs実験群4＆実験群5vs実験群4
mB1_h10_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h10_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h11_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h11_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h12_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h12_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h13_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_self+1),ctl), data=dtmp)
mB1_h13_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_self+1),ctl), data=dtmp)

mB1r_h10_1 <- coeftest(mB1_h10_1, vcov.=vcovHC(mB1_h10_1,"HC2"))
mB1r_h10_2 <- coeftest(mB1_h10_2, vcov.=vcovHC(mB1_h10_2,"HC2"))
mB1r_h11_1 <- coeftest(mB1_h11_1, vcov.=vcovHC(mB1_h11_1,"HC2"))
mB1r_h11_2 <- coeftest(mB1_h11_2, vcov.=vcovHC(mB1_h11_2,"HC2"))
mB1r_h12_1 <- coeftest(mB1_h12_1, vcov.=vcovHC(mB1_h12_1,"HC2"))
mB1r_h12_2 <- coeftest(mB1_h12_2, vcov.=vcovHC(mB1_h12_2,"HC2"))
mB1r_h13_1 <- coeftest(mB1_h13_1, vcov.=vcovHC(mB1_h13_1,"HC2"))
mB1r_h13_2 <- coeftest(mB1_h13_2, vcov.=vcovHC(mB1_h13_2,"HC2"))

mB2_h10_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h10_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h11_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h11_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h12_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h12_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h13_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_self-1),ctl), data=dtmp)
mB2_h13_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_self-1),ctl), data=dtmp)

mB2r_h10_1 <- coeftest(mB2_h10_1, vcov.=vcovHC(mB2_h10_1,"HC2"))
mB2r_h10_2 <- coeftest(mB2_h10_2, vcov.=vcovHC(mB2_h10_2,"HC2"))
mB2r_h11_1 <- coeftest(mB2_h11_1, vcov.=vcovHC(mB2_h11_1,"HC2"))
mB2r_h11_2 <- coeftest(mB2_h11_2, vcov.=vcovHC(mB2_h11_2,"HC2"))
mB2r_h12_1 <- coeftest(mB2_h12_1, vcov.=vcovHC(mB2_h12_1,"HC2"))
mB2r_h12_2 <- coeftest(mB2_h12_2, vcov.=vcovHC(mB2_h12_2,"HC2"))
mB2r_h13_1 <- coeftest(mB2_h13_1, vcov.=vcovHC(mB2_h13_1,"HC2"))
mB2r_h13_2 <- coeftest(mB2_h13_2, vcov.=vcovHC(mB2_h13_2,"HC2"))

# 仮説２の検証（３種類）
# 統制群vs実験群1＆実験群2vs実験群4 ＆ 実験群3vs実験群5
mB1_h20_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h20_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h21_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h21_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h22_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_self+1),ctl), data=dtmp)
mB1_h22_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_self+1),ctl), data=dtmp)

mB1r_h20_1 <- coeftest(mB1_h20_1, vcov.=vcovHC(mB1_h20_1,"HC2"))
mB1r_h20_2 <- coeftest(mB1_h20_2, vcov.=vcovHC(mB1_h20_2,"HC2"))
mB1r_h21_1 <- coeftest(mB1_h21_1, vcov.=vcovHC(mB1_h21_1,"HC2"))
mB1r_h21_2 <- coeftest(mB1_h21_2, vcov.=vcovHC(mB1_h21_2,"HC2"))
mB1r_h22_1 <- coeftest(mB1_h22_1, vcov.=vcovHC(mB1_h22_1,"HC2"))
mB1r_h22_2 <- coeftest(mB1_h22_2, vcov.=vcovHC(mB1_h22_2,"HC2"))

mB2_h20_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h20_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h21_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h21_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h22_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_self-1),ctl), data=dtmp)
mB2_h22_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_self-1),ctl), data=dtmp)

mB2r_h20_1 <- coeftest(mB2_h20_1, vcov.=vcovHC(mB2_h20_1,"HC2"))
mB2r_h20_2 <- coeftest(mB2_h20_2, vcov.=vcovHC(mB2_h20_2,"HC2"))
mB2r_h21_1 <- coeftest(mB2_h21_1, vcov.=vcovHC(mB2_h21_1,"HC2"))
mB2r_h21_2 <- coeftest(mB2_h21_2, vcov.=vcovHC(mB2_h21_2,"HC2"))
mB2r_h22_1 <- coeftest(mB2_h22_1, vcov.=vcovHC(mB2_h22_1,"HC2"))
mB2r_h22_2 <- coeftest(mB2_h22_2, vcov.=vcovHC(mB2_h22_2,"HC2"))

#'
#' ### 自己申告イデオロギーに条件付けされた実験情報刺激の限界効果を用いた仮説検証（図A5）
#'

htest <- data.frame(int = rep(c("左派 (10%=-1)", "右派 (90%=1)"), each=14), 
                    dv = rep(c("生活必需品","その他すべて"), each=7),
                    h = rep(c("H1B","H1B",
                              "H1B","H1B",
                              "H2B","H2B",
                              "H2B"),4),
                    cp = rep(c("2.普遍 - 0.統制",
                               "2.普遍 - 3.選別",
                               "4.逆進+普遍 - 1.逆進",
                               "4.逆進+普遍 - 5.逆進+選別",
                               "1.逆進 - 0.統制",
                               "4.普遍+逆進 - 2.普遍",
                               "5.選別+逆進 - 3.選別"),4),
                    rbind(mB1r_h10_1[2,],mB1r_h11_1[2,],mB1r_h12_1[2,],mB1r_h13_1[2,],
                          mB1r_h20_1[2,],mB1r_h21_1[2,],mB1r_h22_1[2,],
                          mB1r_h10_2[2,],mB1r_h11_2[2,],mB1r_h12_2[2,],mB1r_h13_2[2,],
                          mB1r_h20_2[2,],mB1r_h21_2[2,],mB1r_h22_2[2,],
                          mB2r_h10_1[2,],mB2r_h11_1[2,],mB2r_h12_1[2,],mB2r_h13_1[2,],
                          mB2r_h20_1[2,],mB2r_h21_1[2,],mB2r_h22_1[2,],
                          mB2r_h10_2[2,],mB2r_h11_2[2,],mB2r_h12_2[2,],mB2r_h13_2[2,],
                          mB2r_h20_2[2,],mB2r_h21_2[2,],mB2r_h22_2[2,]))
htest$dv <- factor(htest$dv, levels=unique(htest$dv))
htest$cp <- factor(htest$cp, levels=rev(unique(htest$cp)))
htest$int <- factor(htest$int, levels=unique(htest$int))
levels(htest$int)
htest$lo95 <- htest$Estimate - qnorm(0.975)*htest$Std..Error
htest$up95 <- htest$Estimate + qnorm(0.975)*htest$Std..Error
htest$lo90 <- htest$Estimate - qnorm(0.95)*htest$Std..Error
htest$up90 <- htest$Estimate + qnorm(0.95)*htest$Std..Error

p <- ggplot(htest, aes(x=cp, y=Estimate)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lo95,ymax=up95, color=int), width=0.25, position = position_dodge(width=-0.5)) + 
  geom_errorbar(aes(ymin=lo90,ymax=up90, color=int), width=0, size=1.5, position = position_dodge(width=-0.5)) + 
  geom_point(aes(color=int, shape=int), size=3, position = position_dodge(width=-0.5)) +
  facet_grid(h~dv, scales = "free_y", space = "free_y") + 
  coord_flip() + 
  scale_color_brewer(name="自己申告イデオロギー", type="qual", palette=2) + 
  scale_shape_discrete(name="自己申告イデオロギー") + 
  labs(x=NULL, y="実験刺激効果（従属変数は理想消費税率）", 
       caption="分析の詳細は回帰表を参照。太線は90%信頼区間、細線は95%信頼区間を示している。",
       subtitle = "対象商品") + 
  theme_bw() + theme(legend.position="bottom",
                     plot.subtitle = element_text(hjust=0.5),
                     strip.placement = "outside")
#+ warning=FALSE, fig.width=8, fig.height=5, dev="png", dpi=300, echo=FALSE
p

#+ warning=FALSE, eval=FALSE
# ggsave("htest_mB_originalscale.png", p, width=8, height=6)

#'
#' ## 争点態度イデオロギー条件付け
#'

#'
#' ### 外交安全保障
#' 

mx_ctax1 <- lm(update(tax1_opi ~ as.factor(g_ctax_N)*ide_iss_1,ctl), data=dtmp)
coeftest(mx_ctax1, vcovHC(mx_ctax1, "HC2"))
mBa_1 <- mx_ctax1
mx_ctax2 <- lm(update(tax2_opi ~ as.factor(g_ctax_N)*ide_iss_1,ctl), data=dtmp)
coeftest(mx_ctax2, vcovHC(mx_ctax2, "HC2"))
mBa_2 <- mx_ctax2

# 仮説１の検証（４種類）
# 統制群vs実験群2＆実験群3vs実験群2＆実験群1vs実験群4＆実験群5vs実験群4
mBa1_h10_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h10_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h11_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h11_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h12_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h12_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h13_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h13_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_iss_1+1.35),ctl), data=dtmp)

mBa1r_h10_1 <- coeftest(mBa1_h10_1, vcov.=vcovHC(mBa1_h10_1,"HC2"))
mBa1r_h10_2 <- coeftest(mBa1_h10_2, vcov.=vcovHC(mBa1_h10_2,"HC2"))
mBa1r_h11_1 <- coeftest(mBa1_h11_1, vcov.=vcovHC(mBa1_h11_1,"HC2"))
mBa1r_h11_2 <- coeftest(mBa1_h11_2, vcov.=vcovHC(mBa1_h11_2,"HC2"))
mBa1r_h12_1 <- coeftest(mBa1_h12_1, vcov.=vcovHC(mBa1_h12_1,"HC2"))
mBa1r_h12_2 <- coeftest(mBa1_h12_2, vcov.=vcovHC(mBa1_h12_2,"HC2"))
mBa1r_h13_1 <- coeftest(mBa1_h13_1, vcov.=vcovHC(mBa1_h13_1,"HC2"))
mBa1r_h13_2 <- coeftest(mBa1_h13_2, vcov.=vcovHC(mBa1_h13_2,"HC2"))

mBa2_h10_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h10_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h11_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h11_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h12_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h12_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h13_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h13_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_iss_1-1.53),ctl), data=dtmp)

mBa2r_h10_1 <- coeftest(mBa2_h10_1, vcov.=vcovHC(mBa2_h10_1,"HC2"))
mBa2r_h10_2 <- coeftest(mBa2_h10_2, vcov.=vcovHC(mBa2_h10_2,"HC2"))
mBa2r_h11_1 <- coeftest(mBa2_h11_1, vcov.=vcovHC(mBa2_h11_1,"HC2"))
mBa2r_h11_2 <- coeftest(mBa2_h11_2, vcov.=vcovHC(mBa2_h11_2,"HC2"))
mBa2r_h12_1 <- coeftest(mBa2_h12_1, vcov.=vcovHC(mBa2_h12_1,"HC2"))
mBa2r_h12_2 <- coeftest(mBa2_h12_2, vcov.=vcovHC(mBa2_h12_2,"HC2"))
mBa2r_h13_1 <- coeftest(mBa2_h13_1, vcov.=vcovHC(mBa2_h13_1,"HC2"))
mBa2r_h13_2 <- coeftest(mBa2_h13_2, vcov.=vcovHC(mBa2_h13_2,"HC2"))

# 仮説２の検証（３種類）
# 統制群vs実験群1＆実験群2vs実験群4 ＆ 実験群3vs実験群5
mBa1_h20_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h20_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h21_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h21_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h22_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h22_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_iss_1+1.35),ctl), data=dtmp)

mBa1r_h20_1 <- coeftest(mBa1_h20_1, vcov.=vcovHC(mBa1_h20_1,"HC2"))
mBa1r_h20_2 <- coeftest(mBa1_h20_2, vcov.=vcovHC(mBa1_h20_2,"HC2"))
mBa1r_h21_1 <- coeftest(mBa1_h21_1, vcov.=vcovHC(mBa1_h21_1,"HC2"))
mBa1r_h21_2 <- coeftest(mBa1_h21_2, vcov.=vcovHC(mBa1_h21_2,"HC2"))
mBa1r_h22_1 <- coeftest(mBa1_h22_1, vcov.=vcovHC(mBa1_h22_1,"HC2"))
mBa1r_h22_2 <- coeftest(mBa1_h22_2, vcov.=vcovHC(mBa1_h22_2,"HC2"))

mBa2_h20_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h20_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h21_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h21_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h22_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h22_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_iss_1-1.53),ctl), data=dtmp)

mBa2r_h20_1 <- coeftest(mBa2_h20_1, vcov.=vcovHC(mBa2_h20_1,"HC2"))
mBa2r_h20_2 <- coeftest(mBa2_h20_2, vcov.=vcovHC(mBa2_h20_2,"HC2"))
mBa2r_h21_1 <- coeftest(mBa2_h21_1, vcov.=vcovHC(mBa2_h21_1,"HC2"))
mBa2r_h21_2 <- coeftest(mBa2_h21_2, vcov.=vcovHC(mBa2_h21_2,"HC2"))
mBa2r_h22_1 <- coeftest(mBa2_h22_1, vcov.=vcovHC(mBa2_h22_1,"HC2"))
mBa2r_h22_2 <- coeftest(mBa2_h22_2, vcov.=vcovHC(mBa2_h22_2,"HC2"))

#' 
#' ### 権利機会平等
#' 

mx_ctax1 <- lm(update(tax1_opi ~ as.factor(g_ctax_N)*ide_iss_2,ctl), data=dtmp)
coeftest(mx_ctax1, vcovHC(mx_ctax1, "HC2"))
mBb_1 <- mx_ctax1
mx_ctax2 <- lm(update(tax2_opi ~ as.factor(g_ctax_N)*ide_iss_2,ctl), data=dtmp)
coeftest(mx_ctax2, vcovHC(mx_ctax2, "HC2"))
mBb_2 <- mx_ctax2

# 仮説１の検証（４種類）
# 統制群vs実験群2＆実験群3vs実験群2＆実験群1vs実験群4＆実験群5vs実験群4
mBb1_h10_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h10_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h11_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h11_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h12_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h12_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h13_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h13_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_iss_2+1.50),ctl), data=dtmp)

mBb1r_h10_1 <- coeftest(mBb1_h10_1, vcov.=vcovHC(mBb1_h10_1,"HC2"))
mBb1r_h10_2 <- coeftest(mBb1_h10_2, vcov.=vcovHC(mBb1_h10_2,"HC2"))
mBb1r_h11_1 <- coeftest(mBb1_h11_1, vcov.=vcovHC(mBb1_h11_1,"HC2"))
mBb1r_h11_2 <- coeftest(mBb1_h11_2, vcov.=vcovHC(mBb1_h11_2,"HC2"))
mBb1r_h12_1 <- coeftest(mBb1_h12_1, vcov.=vcovHC(mBb1_h12_1,"HC2"))
mBb1r_h12_2 <- coeftest(mBb1_h12_2, vcov.=vcovHC(mBb1_h12_2,"HC2"))
mBb1r_h13_1 <- coeftest(mBb1_h13_1, vcov.=vcovHC(mBb1_h13_1,"HC2"))
mBb1r_h13_2 <- coeftest(mBb1_h13_2, vcov.=vcovHC(mBb1_h13_2,"HC2"))

mBb2_h10_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h10_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h11_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h11_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h12_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h12_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h13_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h13_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_iss_2-1.48),ctl), data=dtmp)

mBb2r_h10_1 <- coeftest(mBb2_h10_1, vcov.=vcovHC(mBb2_h10_1,"HC2"))
mBb2r_h10_2 <- coeftest(mBb2_h10_2, vcov.=vcovHC(mBb2_h10_2,"HC2"))
mBb2r_h11_1 <- coeftest(mBb2_h11_1, vcov.=vcovHC(mBb2_h11_1,"HC2"))
mBb2r_h11_2 <- coeftest(mBb2_h11_2, vcov.=vcovHC(mBb2_h11_2,"HC2"))
mBb2r_h12_1 <- coeftest(mBb2_h12_1, vcov.=vcovHC(mBb2_h12_1,"HC2"))
mBb2r_h12_2 <- coeftest(mBb2_h12_2, vcov.=vcovHC(mBb2_h12_2,"HC2"))
mBb2r_h13_1 <- coeftest(mBb2_h13_1, vcov.=vcovHC(mBb2_h13_1,"HC2"))
mBb2r_h13_2 <- coeftest(mBb2_h13_2, vcov.=vcovHC(mBb2_h13_2,"HC2"))

# 仮説２の検証（３種類）
# 統制群vs実験群1＆実験群2vs実験群4 ＆ 実験群3vs実験群5
mBb1_h20_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h20_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h21_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h21_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h22_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h22_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_iss_2+1.50),ctl), data=dtmp)

mBb1r_h20_1 <- coeftest(mBb1_h20_1, vcov.=vcovHC(mBb1_h20_1,"HC2"))
mBb1r_h20_2 <- coeftest(mBb1_h20_2, vcov.=vcovHC(mBb1_h20_2,"HC2"))
mBb1r_h21_1 <- coeftest(mBb1_h21_1, vcov.=vcovHC(mBb1_h21_1,"HC2"))
mBb1r_h21_2 <- coeftest(mBb1_h21_2, vcov.=vcovHC(mBb1_h21_2,"HC2"))
mBb1r_h22_1 <- coeftest(mBb1_h22_1, vcov.=vcovHC(mBb1_h22_1,"HC2"))
mBb1r_h22_2 <- coeftest(mBb1_h22_2, vcov.=vcovHC(mBb1_h22_2,"HC2"))

mBb2_h20_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h20_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h21_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h21_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h22_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h22_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_iss_2-1.48),ctl), data=dtmp)

mBb2r_h20_1 <- coeftest(mBb2_h20_1, vcov.=vcovHC(mBb2_h20_1,"HC2"))
mBb2r_h20_2 <- coeftest(mBb2_h20_2, vcov.=vcovHC(mBb2_h20_2,"HC2"))
mBb2r_h21_1 <- coeftest(mBb2_h21_1, vcov.=vcovHC(mBb2_h21_1,"HC2"))
mBb2r_h21_2 <- coeftest(mBb2_h21_2, vcov.=vcovHC(mBb2_h21_2,"HC2"))
mBb2r_h22_1 <- coeftest(mBb2_h22_1, vcov.=vcovHC(mBb2_h22_1,"HC2"))
mBb2r_h22_2 <- coeftest(mBb2_h22_2, vcov.=vcovHC(mBb2_h22_2,"HC2"))

#'
#' ## 政党支持イデオロギー条件付け
#'

mx_ctax1 <- lm(update(tax1_opi ~ as.factor(g_ctax_N)*ide_psup,ctl), data=dtmp)
coeftest(mx_ctax1, vcovHC(mx_ctax1, "HC2"))
mBc_1 <- mx_ctax1
mx_ctax2 <- lm(update(tax2_opi ~ as.factor(g_ctax_N)*ide_psup,ctl), data=dtmp)
coeftest(mx_ctax2, vcovHC(mx_ctax2, "HC2"))
mBc_2 <- mx_ctax2

# 仮説１の検証（４種類）
# 統制群vs実験群2＆実験群3vs実験群2＆実験群1vs実験群4＆実験群5vs実験群4
mBc1_h10_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h10_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h11_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h11_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h12_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h12_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h13_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h13_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_psup+1),ctl), data=dtmp)

mBc1r_h10_1 <- coeftest(mBc1_h10_1, vcov.=vcovHC(mBc1_h10_1,"HC2"))
mBc1r_h10_2 <- coeftest(mBc1_h10_2, vcov.=vcovHC(mBc1_h10_2,"HC2"))
mBc1r_h11_1 <- coeftest(mBc1_h11_1, vcov.=vcovHC(mBc1_h11_1,"HC2"))
mBc1r_h11_2 <- coeftest(mBc1_h11_2, vcov.=vcovHC(mBc1_h11_2,"HC2"))
mBc1r_h12_1 <- coeftest(mBc1_h12_1, vcov.=vcovHC(mBc1_h12_1,"HC2"))
mBc1r_h12_2 <- coeftest(mBc1_h12_2, vcov.=vcovHC(mBc1_h12_2,"HC2"))
mBc1r_h13_1 <- coeftest(mBc1_h13_1, vcov.=vcovHC(mBc1_h13_1,"HC2"))
mBc1r_h13_2 <- coeftest(mBc1_h13_2, vcov.=vcovHC(mBc1_h13_2,"HC2"))

mBc2_h10_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h10_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h11_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h11_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h12_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h12_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h13_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h13_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_psup-1),ctl), data=dtmp)

mBc2r_h10_1 <- coeftest(mBc2_h10_1, vcov.=vcovHC(mBc2_h10_1,"HC2"))
mBc2r_h10_2 <- coeftest(mBc2_h10_2, vcov.=vcovHC(mBc2_h10_2,"HC2"))
mBc2r_h11_1 <- coeftest(mBc2_h11_1, vcov.=vcovHC(mBc2_h11_1,"HC2"))
mBc2r_h11_2 <- coeftest(mBc2_h11_2, vcov.=vcovHC(mBc2_h11_2,"HC2"))
mBc2r_h12_1 <- coeftest(mBc2_h12_1, vcov.=vcovHC(mBc2_h12_1,"HC2"))
mBc2r_h12_2 <- coeftest(mBc2_h12_2, vcov.=vcovHC(mBc2_h12_2,"HC2"))
mBc2r_h13_1 <- coeftest(mBc2_h13_1, vcov.=vcovHC(mBc2_h13_1,"HC2"))
mBc2r_h13_2 <- coeftest(mBc2_h13_2, vcov.=vcovHC(mBc2_h13_2,"HC2"))

# 仮説２の検証（３種類）
# 統制群vs実験群1＆実験群2vs実験群4 ＆ 実験群3vs実験群5
mBc1_h20_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h20_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h21_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h21_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h22_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h22_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_psup+1),ctl), data=dtmp)

mBc1r_h20_1 <- coeftest(mBc1_h20_1, vcov.=vcovHC(mBc1_h20_1,"HC2"))
mBc1r_h20_2 <- coeftest(mBc1_h20_2, vcov.=vcovHC(mBc1_h20_2,"HC2"))
mBc1r_h21_1 <- coeftest(mBc1_h21_1, vcov.=vcovHC(mBc1_h21_1,"HC2"))
mBc1r_h21_2 <- coeftest(mBc1_h21_2, vcov.=vcovHC(mBc1_h21_2,"HC2"))
mBc1r_h22_1 <- coeftest(mBc1_h22_1, vcov.=vcovHC(mBc1_h22_1,"HC2"))
mBc1r_h22_2 <- coeftest(mBc1_h22_2, vcov.=vcovHC(mBc1_h22_2,"HC2"))

mBc2_h20_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h20_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h21_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h21_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h22_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h22_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_psup-1),ctl), data=dtmp)

mBc2r_h20_1 <- coeftest(mBc2_h20_1, vcov.=vcovHC(mBc2_h20_1,"HC2"))
mBc2r_h20_2 <- coeftest(mBc2_h20_2, vcov.=vcovHC(mBc2_h20_2,"HC2"))
mBc2r_h21_1 <- coeftest(mBc2_h21_1, vcov.=vcovHC(mBc2_h21_1,"HC2"))
mBc2r_h21_2 <- coeftest(mBc2_h21_2, vcov.=vcovHC(mBc2_h21_2,"HC2"))
mBc2r_h22_1 <- coeftest(mBc2_h22_1, vcov.=vcovHC(mBc2_h22_1,"HC2"))
mBc2r_h22_2 <- coeftest(mBc2_h22_2, vcov.=vcovHC(mBc2_h22_2,"HC2"))

#'
#' ## 交差項による仮説検証
#' 

htest0 <- data.frame(int = rep(c("世帯収入", "自己申告イデオロギー",
                                 "外交安全保障イデオロギー",
                                 "権利機会平等イデオロギー",
                                 "政党支持イデオロギー"
), each=14),
dv = rep(c("生活必需品","その他すべて"), each=7),
h = rep(c("H1A/B","H1A/B",
          "H1A/B","H1A/B",
          "H2A/B","H2A/B",
          "H2A/B"),2*5),
cp = rep(c("2.普遍 - 0.統制",
           "2.普遍 - 3.選別",
           "4.逆進+普遍 - 1.逆進",
           "4.逆進+普遍 - 5.逆進+選別",
           "1.逆進 - 0.統制",
           "4.普遍+逆進 - 2.普遍",
           "5.選別+逆進 - 3.選別"),2*5),
rbind(mA1r_h10_1[18,],mA1r_h11_1[18,],mA1r_h12_1[18,],mA1r_h13_1[18,],
      mA1r_h20_1[18,],mA1r_h21_1[18,],mA1r_h22_1[18,],
      mA1r_h10_2[18,],mA1r_h11_2[18,],mA1r_h12_2[18,],mA1r_h13_2[18,],
      mA1r_h20_2[18,],mA1r_h21_2[18,],mA1r_h22_2[18,],
      mB1r_h10_1[18,],mB1r_h11_1[18,],mB1r_h12_1[18,],mB1r_h13_1[18,],
      mB1r_h20_1[18,],mB1r_h21_1[18,],mB1r_h22_1[18,],
      mB1r_h10_2[18,],mB1r_h11_2[18,],mB1r_h12_2[18,],mB1r_h13_2[18,],
      mB1r_h20_2[18,],mB1r_h21_2[18,],mB1r_h22_2[18,],
      mBa1r_h10_1[18,],mBa1r_h11_1[18,],mBa1r_h12_1[18,],mBa1r_h13_1[18,],
      mBa1r_h20_1[18,],mBa1r_h21_1[18,],mBa1r_h22_1[18,],
      mBa1r_h10_2[18,],mBa1r_h11_2[18,],mBa1r_h12_2[18,],mBa1r_h13_2[18,],
      mBa1r_h20_2[18,],mBa1r_h21_2[18,],mBa1r_h22_2[18,],
      mBb1r_h10_1[18,],mBb1r_h11_1[18,],mBb1r_h12_1[18,],mBb1r_h13_1[18,],
      mBb1r_h20_1[18,],mBb1r_h21_1[18,],mBb1r_h22_1[18,],
      mBb1r_h10_2[18,],mBb1r_h11_2[18,],mBb1r_h12_2[18,],mBb1r_h13_2[18,],
      mBb1r_h20_2[18,],mBb1r_h21_2[18,],mBb1r_h22_2[18,],
      mBc1r_h10_1[18,],mBc1r_h11_1[18,],mBc1r_h12_1[18,],mBc1r_h13_1[18,],
      mBc1r_h20_1[18,],mBc1r_h21_1[18,],mBc1r_h22_1[18,],
      mBc1r_h10_2[18,],mBc1r_h11_2[18,],mBc1r_h12_2[18,],mBc1r_h13_2[18,],
      mBc1r_h20_2[18,],mBc1r_h21_2[18,],mBc1r_h22_2[18,]))
htest0$int <- factor(htest0$int, levels=unique(htest0$int))
htest0$dv <- factor(htest0$dv, levels=unique(htest0$dv))
htest0$cp <- factor(htest0$cp, levels=rev(unique(htest0$cp)))
htest0$lo95 <- htest0$Estimate - qnorm(0.975)*htest0$Std..Error
htest0$up95 <- htest0$Estimate + qnorm(0.975)*htest0$Std..Error
htest0$lo90 <- htest0$Estimate - qnorm(0.95)*htest0$Std..Error
htest0$up90 <- htest0$Estimate + qnorm(0.95)*htest0$Std..Error

#' 
#' ### 実験情報刺激効果と収入・イデオロギーの交差項係数による仮説検証（図1）
#'   

p <- ggplot(subset(htest0, int%in%c("世帯収入","自己申告イデオロギー")), 
            aes(x=cp, y=Estimate)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lo95,ymax=up95, color=int, alpha=int), 
                width=0.25, position = position_dodge(width=-0.7)) + 
  geom_errorbar(aes(ymin=lo90,ymax=up90, color=int, alpha=int), 
                width=0, size=1.5, position = position_dodge(width=-0.7)) + 
  geom_point(aes(color=int, shape=int, alpha=int), size=3, 
             position = position_dodge(width=-0.7)) +
  facet_grid(h~dv, scales = "free_y", space = "free_y") + 
  coord_flip() + 
  scale_color_manual(name="収入・イデオロギー", values=rep("black",2)) + 
  # scale_color_brewer(name="収入・イデオロギー", type="qual", palette=2) + 
  scale_shape_discrete(name="収入・イデオロギー") + 
  scale_alpha_manual(name = "収入・イデオロギー", values=c(1,rep(0.3,1))) + 
  labs(x=NULL, y="実験群比較変数と収入・イデオロギーの交差項係数\n（従属変数は理想消費税率）", 
       caption="分析の詳細は回帰表を参照。太線は90%信頼区間、細線は95%信頼区間を示している。",
       subtitle = "対象商品") + 
  theme_bw() + theme(legend.position="bottom",
                     plot.subtitle = element_text(hjust=0.5),
                     strip.placement = "outside")

#+ warning=FALSE, fig.width=8, fig.height=6, dev="png", dpi=300, echo=FALSE
p

#+ warning=FALSE, eval=FALSE
# ggsave("htest0_v1_originalscale.png", p, width=8, height=6)

#' 
#' ### 実験情報刺激効果と収入・イデオロギーの交差項係数による仮説検証（その他のイデオロギー指標を含む）（図A8）
#'   

p <- ggplot(htest0, 
            aes(x=cp, y=Estimate)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lo95,ymax=up95, color=int, alpha=int), 
                width=0.25, position = position_dodge(width=-0.7)) + 
  geom_errorbar(aes(ymin=lo90,ymax=up90, color=int, alpha=int), 
                width=0, size=1.5, position = position_dodge(width=-0.7)) + 
  geom_point(aes(color=int, shape=int, alpha=int), size=3, 
             position = position_dodge(width=-0.7)) +
  facet_grid(h~dv, scales = "free_y", space = "free_y") + 
  coord_flip() + 
  scale_color_manual(name="収入・イデオロギー", values=rep("black",5)) + 
  # scale_color_brewer(name="収入・イデオロギー", type="qual", palette=2) + 
  scale_shape_discrete(name="収入・イデオロギー") + 
  scale_alpha_manual(name = "収入・イデオロギー", values=c(1,rep(0.3,4))) + 
  labs(x=NULL, y="実験群比較変数と収入・イデオロギーの交差項係数\n（従属変数は理想消費税率）", 
       caption="分析の詳細は回帰表を参照。太線は90%信頼区間、細線は95%信頼区間を示している。",
       subtitle = "対象商品") + 
  theme_bw() + theme(legend.position="bottom",
                     plot.subtitle = element_text(hjust=0.5),
                     strip.placement = "outside") + 
  guides(color=guide_legend(nrow=2,byrow=TRUE),
         shape=guide_legend(nrow=2,byrow=TRUE),
         alpha=guide_legend(nrow=2,byrow=TRUE))

#+ warning=FALSE, fig.width=9, fig.height=7, dev="png", dpi=300, echo=FALSE
p

#+ warning=FALSE, eval=FALSE
# ggsave("htest0_v2_originalscale.png", p, width=9, height=7)

#'
#' ## 表のエクスポート
#'

vnxinc_rev <- vnxinc[c(1:7,8:17,18:22)]
vnxinc_rev <- gsub("世帯収入","収入／イデオロギー",vnxinc_rev)
vnmap <- list()
for(i in c(1:7,18:22)) vnmap[[names(coef(mA_1))[i]]] = vnxinc_rev[i]
for(i in c(7,18:22)) vnmap[[names(coef(mB_1))[i]]] = vnxinc_rev[i]
for(i in c(7,18:22)) vnmap[[names(coef(mBa_1))[i]]] = vnxinc_rev[i]
for(i in c(7,18:22)) vnmap[[names(coef(mBb_1))[i]]] = vnxinc_rev[i]
for(i in c(7,18:22)) vnmap[[names(coef(mBc_1))[i]]] = vnxinc_rev[i]
for(i in c(8:17)) vnmap[[names(coef(mA_1))[i]]] = vnxinc_rev[i]

vnxA <- gsub("収入／イデオロギー","世帯収入",vnxinc_rev)
vnxB <- gsub("収入／イデオロギー","イデオロギー",vnxinc_rev)
vnmapS <- list()
for(i in c(1:7,18:22)) vnmapS[[names(coef(mA_1))[i]]] = vnxA[i]
for(i in c(7,18:22)) vnmapS[[names(coef(mB_1))[i]]] = vnxB[i]
for(i in c(7,18:22)) vnmapS[[names(coef(mBa_1))[i]]] = vnxB[i]
for(i in c(7,18:22)) vnmapS[[names(coef(mBb_1))[i]]] = vnxB[i]
for(i in c(7,18:22)) vnmapS[[names(coef(mBc_1))[i]]] = vnxB[i]
for(i in c(8:17)) vnmapS[[names(coef(mA_1))[i]]] = vnxA[i]

vnmap2 <- vnmap
for(k in names(vnmap2)) vnmap2[[k]] <- gsub("収入／イデオロギー", "イデオロギー", vnmap2[[k]])

#'
#' ### 直接効果
#'

## 直接効果
screenreg(list(m0_1,m0_2), 
          override.se = list(coeftest(m0_1,vcovHC(m0_1,"HC2"))[,2],
                             coeftest(m0_2,vcovHC(m0_2,"HC2"))[,2]),
          override.pvalues = list(coeftest(m0_1,vcovHC(m0_1,"HC2"))[,4],
                                  coeftest(m0_2,vcovHC(m0_2,"HC2"))[,4]),
          symbol = "+",
          single.row=TRUE, digits = 3, stars = c(0.001,0.01,0.05,0.1),
          custom.coef.map = vnmap,
          custom.model.names = c("1:生活必需品","2:その他すべて"),
          caption = "理想消費税率に実験情報刺激が与える効果（重回帰分析）",
          caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
          label="basetab_os", dcolumn = TRUE, booktabs = TRUE, use.packages = FALSE,
          custom.note = "%stars. 最小二乗法による重回帰分析、ロバスト標準誤差使用．")
# texreg(list(m0_1,m0_2), 
#        override.se = list(coeftest(m0_1,vcovHC(m0_1,"HC2"))[,2],
#                           coeftest(m0_2,vcovHC(m0_2,"HC2"))[,2]),
#        override.pvalues = list(coeftest(m0_1,vcovHC(m0_1,"HC2"))[,4],
#                                coeftest(m0_2,vcovHC(m0_2,"HC2"))[,4]),
#        # file = "basetab_originalscale.html", symbol = "&dagger;", 
#        file = "basetab_originalscale.tex", symbol = "\\dagger", 
#        single.row=TRUE, digits = 3, stars = c(0.001,0.01,0.05,0.1),
#        custom.coef.map = vnmap,
#        custom.model.names = c("1:生活必需品","2:その他すべて"),
#        caption = "理想消費税率に実験情報刺激が与える効果（重回帰分析）",
#        caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
#        label="basetab_os", dcolumn = TRUE, booktabs = TRUE, use.packages = FALSE,
#        custom.note = "%stars. 最小二乗法による重回帰分析、ロバスト標準誤差使用．")
# tmp <- readLines("basetab_originalscale.tex")
# tmp <- gsub("{dagger}","{\\dagger}", tmp, fixed=TRUE)
# writeLines(tmp, "basetab_originalscale.tex", useBytes = TRUE)

#'
#' ### 収入と自己申告イデオロギーによる条件付け（表1・表A1）
#'

## Simplified Table（表1）
screenreg(list(mA_1,mA_2,mB_1,mB_2), 
          override.se = list(coeftest(mA_1,vcovHC(mA_1,"HC2"))[,2],
                             coeftest(mA_2,vcovHC(mA_2,"HC2"))[,2],
                             coeftest(mB_1,vcovHC(mB_1,"HC2"))[,2],
                             coeftest(mB_2,vcovHC(mB_2,"HC2"))[,2]),
          override.pvalues = list(coeftest(mA_1,vcovHC(mA_1,"HC2"))[,4],
                                  coeftest(mA_2,vcovHC(mA_2,"HC2"))[,4],
                                  coeftest(mB_1,vcovHC(mB_1,"HC2"))[,4],
                                  coeftest(mB_2,vcovHC(mB_2,"HC2"))[,4]),
          symbol = "+", 
          single.row=FALSE, digits = 3, stars = c(0.001,0.01,0.05,0.1),
          custom.coef.map = vnmapS[c(2:18)],
          custom.model.names = c("1:生活必需品","2:その他すべて",
                                 "3:生活必需品","4:その他すべて"),
          custom.header = list("世帯収入" = 1:2, "自己申告イデオロギー" = 3:4), 
          caption = "理想消費税率に実験情報刺激が与える効果と収入・イデオロギー（重回帰分析）",
          caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
          label="maintab_os", dcolumn = TRUE, booktabs = TRUE, use.packages = FALSE,
          custom.note = "%stars. 最小二乗法による重回帰分析、ロバスト標準誤差使用．統制変数の係数はオンライン補遺を参照．")
# texreg(list(mA_1,mA_2,mB_1,mB_2), 
#        override.se = list(coeftest(mA_1,vcovHC(mA_1,"HC2"))[,2],
#                           coeftest(mA_2,vcovHC(mA_2,"HC2"))[,2],
#                           coeftest(mB_1,vcovHC(mB_1,"HC2"))[,2],
#                           coeftest(mB_2,vcovHC(mB_2,"HC2"))[,2]),
#        override.pvalues = list(coeftest(mA_1,vcovHC(mA_1,"HC2"))[,4],
#                                coeftest(mA_2,vcovHC(mA_2,"HC2"))[,4],
#                                coeftest(mB_1,vcovHC(mB_1,"HC2"))[,4],
#                                coeftest(mB_2,vcovHC(mB_2,"HC2"))[,4]),
#        # file = "maintab_originalscale.html", symbol = "&dagger;",
#        file = "maintab_originalscale.tex", symbol = "\\dagger", 
#        single.row=FALSE, digits = 3, stars = c(0.001,0.01,0.05,0.1),
#        custom.coef.map = vnmapS[c(2:18)],
#        custom.model.names = c("1:生活必需品","2:その他すべて",
#                               "3:生活必需品","4:その他すべて"),
#        custom.header = list("世帯収入" = 1:2, "自己申告イデオロギー" = 3:4), 
#        caption = "理想消費税率に実験情報刺激が与える効果と収入・イデオロギー（重回帰分析）",
#        caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
#        label="maintab_os", dcolumn = TRUE, booktabs = TRUE, use.packages = FALSE,
#        custom.note = "%stars. 最小二乗法による重回帰分析、ロバスト標準誤差使用．統制変数の係数はオンライン補遺を参照．")
# tmp <- readLines("maintab_originalscale.tex")
# tmp <- gsub("{dagger}","{\\dagger}", tmp, fixed=TRUE)
# writeLines(tmp, "maintab_originalscale.tex", useBytes = TRUE)


## Full Table（表A1）
screenreg(list(mA_1,mA_2,mB_1,mB_2), 
          override.se = list(coeftest(mA_1,vcovHC(mA_1,"HC2"))[,2],
                             coeftest(mA_2,vcovHC(mA_2,"HC2"))[,2],
                             coeftest(mB_1,vcovHC(mB_1,"HC2"))[,2],
                             coeftest(mB_2,vcovHC(mB_2,"HC2"))[,2]),
          override.pvalues = list(coeftest(mA_1,vcovHC(mA_1,"HC2"))[,4],
                                  coeftest(mA_2,vcovHC(mA_2,"HC2"))[,4],
                                  coeftest(mB_1,vcovHC(mB_1,"HC2"))[,4],
                                  coeftest(mB_2,vcovHC(mB_2,"HC2"))[,4]),
          symbol = "+", 
          single.row=FALSE, digits = 3, stars = c(0.001,0.01,0.05,0.1),
          custom.coef.map = vnmap,
          custom.model.names = c("1:生活必需品","2:その他すべて",
                                 "3:生活必需品","4:その他すべて"),
          custom.header = list("世帯収入" = 1:2, "自己申告イデオロギー" = 3:4), 
          caption = "理想消費税率に実験情報刺激が与える効果と収入・イデオロギー（重回帰分析）",
          caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
          label="maintab_full_os", dcolumn = TRUE, booktabs = TRUE, use.packages = FALSE,
          custom.note = "%stars. 最小二乗法による重回帰分析、ロバスト標準誤差使用．")
# texreg(list(mA_1,mA_2,mB_1,mB_2), 
#         override.se = list(coeftest(mA_1,vcovHC(mA_1,"HC2"))[,2],
#                            coeftest(mA_2,vcovHC(mA_2,"HC2"))[,2],
#                            coeftest(mB_1,vcovHC(mB_1,"HC2"))[,2],
#                            coeftest(mB_2,vcovHC(mB_2,"HC2"))[,2]),
#         override.pvalues = list(coeftest(mA_1,vcovHC(mA_1,"HC2"))[,4],
#                                 coeftest(mA_2,vcovHC(mA_2,"HC2"))[,4],
#                                 coeftest(mB_1,vcovHC(mB_1,"HC2"))[,4],
#                                 coeftest(mB_2,vcovHC(mB_2,"HC2"))[,4]),
#         # file = "maintab_full_originalscale.html", symbol = "&dagger;", 
#         file = "maintab_full_originalscale.tex", symbol = "\\dagger", 
#         single.row=FALSE, digits = 3, stars = c(0.001,0.01,0.05,0.1),
#         custom.coef.map = vnmap,
#         custom.model.names = c("1:生活必需品","2:その他すべて",
#                                "3:生活必需品","4:その他すべて"),
#         custom.header = list("世帯収入" = 1:2, "自己申告イデオロギー" = 3:4), 
#         caption = "理想消費税率に実験情報刺激が与える効果と収入・イデオロギー（重回帰分析）",
#         caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
#         label="maintab_full_os", dcolumn = TRUE, booktabs = TRUE, use.packages = FALSE,
#         custom.note = "%stars. 最小二乗法による重回帰分析、ロバスト標準誤差使用．")
# tmp <- readLines("maintab_full_originalscale.tex")
# tmp <- gsub("{dagger}","{\\dagger}", tmp, fixed=TRUE)
# writeLines(tmp, "maintab_full_originalscale.tex", useBytes = TRUE)

#'
#' ### その他のイデオロギーによる条件付け
#'

## 他のイデオロギー（生活必需品）
screenreg(list(mBc_1,mBa_1,mBb_1), 
          override.se = list(coeftest(mBc_1,vcovHC(mBc_1,"HC2"))[,2],
                             coeftest(mBa_1,vcovHC(mBa_1,"HC2"))[,2],
                             coeftest(mBb_1,vcovHC(mBb_1,"HC2"))[,2]),
          override.pvalues = list(coeftest(mBc_1,vcovHC(mBc_1,"HC2"))[,4],
                                  coeftest(mBa_1,vcovHC(mBa_1,"HC2"))[,4],
                                  coeftest(mBb_1,vcovHC(mBb_1,"HC2"))[,4]),
          symbol = "+",
          single.row=TRUE, digits = 3, stars = c(0.001,0.01,0.05,0.1),
          custom.coef.map = vnmap2,
          custom.model.names = c("政党支持","外交安全保障","権利機会平等"),
          caption = "生活必需品の理想消費税率に実験情報刺激が与える効果とイデオロギー（重回帰分析）",
          caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
          label="idetab1_os", dcolumn = TRUE, booktabs = TRUE, use.packages = FALSE,
          custom.note = "%stars. 最小二乗法による重回帰分析、ロバスト標準誤差使用．")
# texreg(list(mBc_1,mBa_1,mBb_1), 
#        override.se = list(coeftest(mBc_1,vcovHC(mBc_1,"HC2"))[,2],
#                           coeftest(mBa_1,vcovHC(mBa_1,"HC2"))[,2],
#                           coeftest(mBb_1,vcovHC(mBb_1,"HC2"))[,2]),
#        override.pvalues = list(coeftest(mBc_1,vcovHC(mBc_1,"HC2"))[,4],
#                                coeftest(mBa_1,vcovHC(mBa_1,"HC2"))[,4],
#                                coeftest(mBb_1,vcovHC(mBb_1,"HC2"))[,4]),
#        # file = "idetab1_originalscale.html", symbol = "&dagger;",
#        file = "idetab1_originalscale.tex", symbol = "\\dagger",
#        single.row=TRUE, digits = 3, stars = c(0.001,0.01,0.05,0.1),
#        custom.coef.map = vnmap2,
#        custom.model.names = c("政党支持","外交安全保障","権利機会平等"),
#        caption = "生活必需品の理想消費税率に実験情報刺激が与える効果とイデオロギー（重回帰分析）",
#        caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
#        label="idetab1_os", dcolumn = TRUE, booktabs = TRUE, use.packages = FALSE,
#        custom.note = "%stars. 最小二乗法による重回帰分析、ロバスト標準誤差使用．")
# tmp <- readLines("idetab1_originalscale.tex")
# tmp <- gsub("{dagger}","{\\dagger}", tmp, fixed=TRUE)
# writeLines(tmp, "idetab1_originalscale.tex", useBytes = TRUE)

## 他のイデオロギー（その他の商品）
screenreg(list(mBc_2,mBa_2,mBb_2), 
          override.se = list(coeftest(mBc_2,vcovHC(mBc_2,"HC2"))[,2],
                             coeftest(mBa_2,vcovHC(mBa_2,"HC2"))[,2],
                             coeftest(mBb_2,vcovHC(mBb_2,"HC2"))[,2]),
          override.pvalues = list(coeftest(mBc_2,vcovHC(mBc_2,"HC2"))[,4],
                                  coeftest(mBa_2,vcovHC(mBa_2,"HC2"))[,4],
                                  coeftest(mBb_2,vcovHC(mBb_2,"HC2"))[,4]),
          symbol = "+",
          single.row=TRUE, digits = 3, stars = c(0.001,0.01,0.05,0.1),
          custom.coef.map = vnmap2,
          custom.model.names = c("政党支持","外交安全保障","権利機会平等"),
          caption = "その他全ての商品の理想消費税率に実験情報刺激が与える効果とイデオロギー（重回帰分析）",
          caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
          label="idetab2_os", dcolumn = TRUE, booktabs = TRUE, use.packages = FALSE,
          custom.note = "%stars. 最小二乗法による重回帰分析、ロバスト標準誤差使用．")
# texreg(list(mBc_2,mBa_2,mBb_2), 
#        override.se = list(coeftest(mBc_2,vcovHC(mBc_2,"HC2"))[,2],
#                           coeftest(mBa_2,vcovHC(mBa_2,"HC2"))[,2],
#                           coeftest(mBb_2,vcovHC(mBb_2,"HC2"))[,2]),
#        override.pvalues = list(coeftest(mBc_2,vcovHC(mBc_2,"HC2"))[,4],
#                                coeftest(mBa_2,vcovHC(mBa_2,"HC2"))[,4],
#                                coeftest(mBb_2,vcovHC(mBb_2,"HC2"))[,4]),
#        # file = "idetab2_originalscale.html", symbol = "&dagger;",
#        file = "idetab2_originalscale.tex", symbol = "\\dagger",
#        single.row=TRUE, digits = 3, stars = c(0.001,0.01,0.05,0.1),
#        custom.coef.map = vnmap2,
#        custom.model.names = c("政党支持","外交安全保障","権利機会平等"),
#        caption = "その他全ての商品の理想消費税率に実験情報刺激が与える効果とイデオロギー（重回帰分析）",
#        caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
#        label="idetab2_os", dcolumn = TRUE, booktabs = TRUE, use.packages = FALSE,
#        custom.note = "%stars. 最小二乗法による重回帰分析、ロバスト標準誤差使用．")
# tmp <- readLines("idetab2_originalscale.tex")
# tmp <- gsub("{dagger}","{\\dagger}", tmp, fixed=TRUE)
# writeLines(tmp, "idetab2_originalscale.tex", useBytes = TRUE)

#'
#' # 実験群比較 (統制変数に政治知識なし)
#'
#' ## 準備

# 統制変数
ctl <- formula( ~ .+ fem + age + lvlen + ownh + 
                  as.factor(edu3) + wk + mar + cld)

#'
#' ## 実験刺激の直接効果
#'

m_ctax1 <- lm(update(tax1_opi ~ as.factor(g_ctax_N),ctl), data=dtmp)
m0_1 <- m_ctax1
m_ctax2 <- lm(update(tax2_opi ~ as.factor(g_ctax_N),ctl), data=dtmp)
m0_2 <- m_ctax2

coeftest(m_ctax1, vcov.=vcovHC(m_ctax1,"HC2"))
coeftest(m_ctax2, vcov.=vcovHC(m_ctax2,"HC2"))



# 普遍性刺激効果の検証（４種類）
# 統制群vs実験群2＆実験群3vs実験群2＆実験群1vs実験群4＆実験群5vs実験群4
m0_h10_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5)),ctl), data=dtmp)
m0_h10_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5)),ctl), data=dtmp)
m0_h11_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5)),ctl), data=dtmp)
m0_h11_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5)),ctl), data=dtmp)
m0_h12_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5)),ctl), data=dtmp)
m0_h12_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5)),ctl), data=dtmp)
m0_h13_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3)),ctl), data=dtmp)
m0_h13_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3)),ctl), data=dtmp)

m0r_h10_1 <- coeftest(m0_h10_1, vcov.=vcovHC(m0_h10_1,"HC2"))
m0r_h10_2 <- coeftest(m0_h10_2, vcov.=vcovHC(m0_h10_2,"HC2"))
m0r_h11_1 <- coeftest(m0_h11_1, vcov.=vcovHC(m0_h11_1,"HC2"))
m0r_h11_2 <- coeftest(m0_h11_2, vcov.=vcovHC(m0_h11_2,"HC2"))
m0r_h12_1 <- coeftest(m0_h12_1, vcov.=vcovHC(m0_h12_1,"HC2"))
m0r_h12_2 <- coeftest(m0_h12_2, vcov.=vcovHC(m0_h12_2,"HC2"))
m0r_h13_1 <- coeftest(m0_h13_1, vcov.=vcovHC(m0_h13_1,"HC2"))
m0r_h13_2 <- coeftest(m0_h13_2, vcov.=vcovHC(m0_h13_2,"HC2"))

# 逆進性刺激効果の検証（３種類）
# 統制群vs実験群1＆実験群2vs実験群4 ＆ 実験群3vs実験群5
m0_h20_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5)),ctl), data=dtmp)
m0_h20_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5)),ctl), data=dtmp)
m0_h21_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5)),ctl), data=dtmp)
m0_h21_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5)),ctl), data=dtmp)
m0_h22_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4)),ctl), data=dtmp)
m0_h22_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4)),ctl), data=dtmp)

m0r_h20_1 <- coeftest(m0_h20_1, vcov.=vcovHC(m0_h20_1,"HC2"))
m0r_h20_2 <- coeftest(m0_h20_2, vcov.=vcovHC(m0_h20_2,"HC2"))
m0r_h21_1 <- coeftest(m0_h21_1, vcov.=vcovHC(m0_h21_1,"HC2"))
m0r_h21_2 <- coeftest(m0_h21_2, vcov.=vcovHC(m0_h21_2,"HC2"))
m0r_h22_1 <- coeftest(m0_h22_1, vcov.=vcovHC(m0_h22_1,"HC2"))
m0r_h22_2 <- coeftest(m0_h22_2, vcov.=vcovHC(m0_h22_2,"HC2"))

htest <- data.frame(dv = rep(c("生活必需品","その他すべて"), each=7),
                    h = rep(c("社会保障普遍性","社会保障普遍性",
                              "社会保障普遍性","社会保障普遍性",
                              "消費税逆進性","消費税逆進性",
                              "消費税逆進性"),2),
                    cp = rep(c("2.普遍 - 0.統制",
                               "2.普遍 - 3.選別",
                               "4.逆進+普遍 - 1.逆進",
                               "4.逆進+普遍 - 5.逆進+選別",
                               "1.逆進 - 0.統制",
                               "4.普遍+逆進 - 2.普遍",
                               "5.選別+逆進 - 3.選別"),2),
                    rbind(m0r_h10_1[2,],m0r_h11_1[2,],m0r_h12_1[2,],m0r_h13_1[2,],
                          m0r_h20_1[2,],m0r_h21_1[2,],m0r_h22_1[2,],
                          m0r_h10_2[2,],m0r_h11_2[2,],m0r_h12_2[2,],m0r_h13_2[2,],
                          m0r_h20_2[2,],m0r_h21_2[2,],m0r_h22_2[2,]))
htest$dv <- factor(htest$dv, levels=unique(htest$dv))
htest$cp <- factor(htest$cp, levels=rev(unique(htest$cp)))
htest$h <- factor(htest$h, levels=unique(htest$h))
htest$lo95 <- htest$Estimate - qnorm(0.975)*htest$Std..Error
htest$up95 <- htest$Estimate + qnorm(0.975)*htest$Std..Error
htest$lo90 <- htest$Estimate - qnorm(0.95)*htest$Std..Error
htest$up90 <- htest$Estimate + qnorm(0.95)*htest$Std..Error

#' 
#' ### 仮説に関係する実験群比較に関連する直接効果（図A15）
#'

p <- ggplot(htest, aes(x=cp, y=Estimate)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lo95,ymax=up95), width=0.25) + 
  geom_errorbar(aes(ymin=lo90,ymax=up90), width=0, size=1.5) + 
  geom_point(size=3) +
  facet_grid(h~dv, scale="free_y", switch="y") + 
  coord_flip() + 
  # scale_color_brewer(name="対象商品", type="qual", palette=2) + 
  labs(x=NULL, y="実験刺激効果（従属変数は理想消費税率）", 
       caption="分析の詳細は回帰表を参照。統制変数有。太線は90%信頼区間、細線は95%信頼区間を示している。", 
       subtitle = "対象商品") + 
  theme_bw() + 
  theme(plot.subtitle = element_text(hjust=0.5),
        strip.background.y = element_blank(),
        strip.text.y = element_blank(),
        strip.placement = "outside")

#+ warning=FALSE, fig.width=8, fig.height=5, dev="png", dpi=300, echo=FALSE
p

#+ warning=FALSE, eval=FALSE
# ggsave("htest_m0_originalscale_wokn.png", p, width=8, height=5)

#'
#' ## 世帯収入条件付け
#'

mx_ctax1 <- lm(update(tax1_opi ~ as.factor(g_ctax_N)*inc,ctl), data=dtmp)
coeftest(mx_ctax1, vcovHC(mx_ctax1, "HC2"))
mA_1 <- mx_ctax1
mx_ctax2 <- lm(update(tax2_opi ~ as.factor(g_ctax_N)*inc,ctl), data=dtmp)
coeftest(mx_ctax2, vcovHC(mx_ctax2, "HC2"))
mA_2 <- mx_ctax2

# 仮説１の検証（４種類）
# 統制群vs実験群2＆実験群3vs実験群2＆実験群1vs実験群4＆実験群5vs実験群4
mA1_h10_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(inc-2),ctl), data=dtmp)
mA1_h10_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(inc-2),ctl), data=dtmp)
mA1_h11_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(inc-2),ctl), data=dtmp)
mA1_h11_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(inc-2),ctl), data=dtmp)
mA1_h12_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(inc-2),ctl), data=dtmp)
mA1_h12_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(inc-2),ctl), data=dtmp)
mA1_h13_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(inc-2),ctl), data=dtmp)
mA1_h13_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(inc-2),ctl), data=dtmp)

mA1r_h10_1 <- coeftest(mA1_h10_1, vcov.=vcovHC(mA1_h10_1,"HC2"))
mA1r_h10_2 <- coeftest(mA1_h10_2, vcov.=vcovHC(mA1_h10_2,"HC2"))
mA1r_h11_1 <- coeftest(mA1_h11_1, vcov.=vcovHC(mA1_h11_1,"HC2"))
mA1r_h11_2 <- coeftest(mA1_h11_2, vcov.=vcovHC(mA1_h11_2,"HC2"))
mA1r_h12_1 <- coeftest(mA1_h12_1, vcov.=vcovHC(mA1_h12_1,"HC2"))
mA1r_h12_2 <- coeftest(mA1_h12_2, vcov.=vcovHC(mA1_h12_2,"HC2"))
mA1r_h13_1 <- coeftest(mA1_h13_1, vcov.=vcovHC(mA1_h13_1,"HC2"))
mA1r_h13_2 <- coeftest(mA1_h13_2, vcov.=vcovHC(mA1_h13_2,"HC2"))

mA2_h10_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(inc-5),ctl), data=dtmp)
mA2_h10_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(inc-5),ctl), data=dtmp)
mA2_h11_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(inc-5),ctl), data=dtmp)
mA2_h11_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(inc-5),ctl), data=dtmp)
mA2_h12_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(inc-5),ctl), data=dtmp)
mA2_h12_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(inc-5),ctl), data=dtmp)
mA2_h13_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(inc-5),ctl), data=dtmp)
mA2_h13_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(inc-5),ctl), data=dtmp)

mA2r_h10_1 <- coeftest(mA2_h10_1, vcov.=vcovHC(mA2_h10_1,"HC2"))
mA2r_h10_2 <- coeftest(mA2_h10_2, vcov.=vcovHC(mA2_h10_2,"HC2"))
mA2r_h11_1 <- coeftest(mA2_h11_1, vcov.=vcovHC(mA2_h11_1,"HC2"))
mA2r_h11_2 <- coeftest(mA2_h11_2, vcov.=vcovHC(mA2_h11_2,"HC2"))
mA2r_h12_1 <- coeftest(mA2_h12_1, vcov.=vcovHC(mA2_h12_1,"HC2"))
mA2r_h12_2 <- coeftest(mA2_h12_2, vcov.=vcovHC(mA2_h12_2,"HC2"))
mA2r_h13_1 <- coeftest(mA2_h13_1, vcov.=vcovHC(mA2_h13_1,"HC2"))
mA2r_h13_2 <- coeftest(mA2_h13_2, vcov.=vcovHC(mA2_h13_2,"HC2"))

# 仮説２の検証（３種類）
# 統制群vs実験群1＆実験群2vs実験群4 ＆ 実験群3vs実験群5
mA1_h20_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(inc-2),ctl), data=dtmp)
mA1_h20_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(inc-2),ctl), data=dtmp)
mA1_h21_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(inc-2),ctl), data=dtmp)
mA1_h21_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(inc-2),ctl), data=dtmp)
mA1_h22_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(inc-2),ctl), data=dtmp)
mA1_h22_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(inc-2),ctl), data=dtmp)

mA1r_h20_1 <- coeftest(mA1_h20_1, vcov.=vcovHC(mA1_h20_1,"HC2"))
mA1r_h20_2 <- coeftest(mA1_h20_2, vcov.=vcovHC(mA1_h20_2,"HC2"))
mA1r_h21_1 <- coeftest(mA1_h21_1, vcov.=vcovHC(mA1_h21_1,"HC2"))
mA1r_h21_2 <- coeftest(mA1_h21_2, vcov.=vcovHC(mA1_h21_2,"HC2"))
mA1r_h22_1 <- coeftest(mA1_h22_1, vcov.=vcovHC(mA1_h22_1,"HC2"))
mA1r_h22_2 <- coeftest(mA1_h22_2, vcov.=vcovHC(mA1_h22_2,"HC2"))

mA2_h20_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(inc-5),ctl), data=dtmp)
mA2_h20_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(inc-5),ctl), data=dtmp)
mA2_h21_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(inc-5),ctl), data=dtmp)
mA2_h21_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(inc-5),ctl), data=dtmp)
mA2_h22_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(inc-5),ctl), data=dtmp)
mA2_h22_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(inc-5),ctl), data=dtmp)

mA2r_h20_1 <- coeftest(mA2_h20_1, vcov.=vcovHC(mA2_h20_1,"HC2"))
mA2r_h20_2 <- coeftest(mA2_h20_2, vcov.=vcovHC(mA2_h20_2,"HC2"))
mA2r_h21_1 <- coeftest(mA2_h21_1, vcov.=vcovHC(mA2_h21_1,"HC2"))
mA2r_h21_2 <- coeftest(mA2_h21_2, vcov.=vcovHC(mA2_h21_2,"HC2"))
mA2r_h22_1 <- coeftest(mA2_h22_1, vcov.=vcovHC(mA2_h22_1,"HC2"))
mA2r_h22_2 <- coeftest(mA2_h22_2, vcov.=vcovHC(mA2_h22_2,"HC2"))

#'
#' ### 世帯収入に条件付けされた実験情報刺激の限界効果を用いた仮説検証（図A14）
#'

htest <- data.frame(int = rep(c("200〜400万円(10%=2)","800〜1000万円(90%=5)"), each=14), 
                    dv = rep(c("生活必需品","その他すべて"), each=7),
                    h = rep(c("H1A","H1A",
                              "H1A","H1A",
                              "H2A","H2A",
                              "H2A"),4),
                    cp = rep(c("2.普遍 - 0.統制",
                               "2.普遍 - 3.選別",
                               "4.逆進+普遍 - 1.逆進",
                               "4.逆進+普遍 - 5.逆進+選別",
                               "1.逆進 - 0.統制",
                               "4.普遍+逆進 - 2.普遍",
                               "5.選別+逆進 - 3.選別"),4),
                    rbind(mA1r_h10_1[2,],mA1r_h11_1[2,],mA1r_h12_1[2,],mA1r_h13_1[2,],
                          mA1r_h20_1[2,],mA1r_h21_1[2,],mA1r_h22_1[2,],
                          mA1r_h10_2[2,],mA1r_h11_2[2,],mA1r_h12_2[2,],mA1r_h13_2[2,],
                          mA1r_h20_2[2,],mA1r_h21_2[2,],mA1r_h22_2[2,],
                          mA2r_h10_1[2,],mA2r_h11_1[2,],mA2r_h12_1[2,],mA2r_h13_1[2,],
                          mA2r_h20_1[2,],mA2r_h21_1[2,],mA2r_h22_1[2,],
                          mA2r_h10_2[2,],mA2r_h11_2[2,],mA2r_h12_2[2,],mA2r_h13_2[2,],
                          mA2r_h20_2[2,],mA2r_h21_2[2,],mA2r_h22_2[2,]))
htest$dv <- factor(htest$dv, levels=unique(htest$dv))
htest$cp <- factor(htest$cp, levels=rev(unique(htest$cp)))
htest$lo95 <- htest$Estimate - qnorm(0.975)*htest$Std..Error
htest$up95 <- htest$Estimate + qnorm(0.975)*htest$Std..Error
htest$lo90 <- htest$Estimate - qnorm(0.95)*htest$Std..Error
htest$up90 <- htest$Estimate + qnorm(0.95)*htest$Std..Error

p <- ggplot(htest, aes(x=cp, y=Estimate)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lo95,ymax=up95, color=int), width=0.25, position = position_dodge(width=-0.5)) + 
  geom_errorbar(aes(ymin=lo90,ymax=up90, color=int), width=0, size=1.5, position = position_dodge(width=-0.5)) + 
  geom_point(aes(color=int, shape=int), size=3, position = position_dodge(width=-0.5)) +
  facet_grid(h~dv, scales = "free_y", space = "free_y") + 
  coord_flip() + 
  scale_color_brewer(name="世帯収入", type="qual", palette=2) + 
  scale_shape_discrete(name="世帯収入") + 
  labs(x=NULL, y="実験刺激効果（従属変数は理想消費税率）", 
       caption="分析の詳細は回帰表を参照。太線は90%信頼区間、細線は95%信頼区間を示している。",
       subtitle = "対象商品") + 
  theme_bw() + theme(legend.position="bottom",
                     plot.subtitle = element_text(hjust=0.5),
                     strip.placement = "outside")
#+ warning=FALSE, fig.width=8, fig.height=5, dev="png", dpi=300, echo=FALSE
p

#+ warning=FALSE, eval=FALSE
# ggsave("htest_mA_originalscale_wokn.png", p, width=8, height=6)

#'
#' ## 自己申告イデオロギー条件付け
#'

mx_ctax1 <- lm(update(tax1_opi ~ as.factor(g_ctax_N)*ide_self,ctl), data=dtmp)
coeftest(mx_ctax1, vcovHC(mx_ctax1, "HC2"))
mB_1 <- mx_ctax1
mx_ctax2 <- lm(update(tax2_opi ~ as.factor(g_ctax_N)*ide_self,ctl), data=dtmp)
coeftest(mx_ctax2, vcovHC(mx_ctax2, "HC2"))
mB_2 <- mx_ctax2

# 仮説１の検証（４種類）
# 統制群vs実験群2＆実験群3vs実験群2＆実験群1vs実験群4＆実験群5vs実験群4
mB1_h10_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h10_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h11_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h11_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h12_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h12_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h13_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_self+1),ctl), data=dtmp)
mB1_h13_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_self+1),ctl), data=dtmp)

mB1r_h10_1 <- coeftest(mB1_h10_1, vcov.=vcovHC(mB1_h10_1,"HC2"))
mB1r_h10_2 <- coeftest(mB1_h10_2, vcov.=vcovHC(mB1_h10_2,"HC2"))
mB1r_h11_1 <- coeftest(mB1_h11_1, vcov.=vcovHC(mB1_h11_1,"HC2"))
mB1r_h11_2 <- coeftest(mB1_h11_2, vcov.=vcovHC(mB1_h11_2,"HC2"))
mB1r_h12_1 <- coeftest(mB1_h12_1, vcov.=vcovHC(mB1_h12_1,"HC2"))
mB1r_h12_2 <- coeftest(mB1_h12_2, vcov.=vcovHC(mB1_h12_2,"HC2"))
mB1r_h13_1 <- coeftest(mB1_h13_1, vcov.=vcovHC(mB1_h13_1,"HC2"))
mB1r_h13_2 <- coeftest(mB1_h13_2, vcov.=vcovHC(mB1_h13_2,"HC2"))

mB2_h10_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h10_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h11_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h11_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h12_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h12_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h13_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_self-1),ctl), data=dtmp)
mB2_h13_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_self-1),ctl), data=dtmp)

mB2r_h10_1 <- coeftest(mB2_h10_1, vcov.=vcovHC(mB2_h10_1,"HC2"))
mB2r_h10_2 <- coeftest(mB2_h10_2, vcov.=vcovHC(mB2_h10_2,"HC2"))
mB2r_h11_1 <- coeftest(mB2_h11_1, vcov.=vcovHC(mB2_h11_1,"HC2"))
mB2r_h11_2 <- coeftest(mB2_h11_2, vcov.=vcovHC(mB2_h11_2,"HC2"))
mB2r_h12_1 <- coeftest(mB2_h12_1, vcov.=vcovHC(mB2_h12_1,"HC2"))
mB2r_h12_2 <- coeftest(mB2_h12_2, vcov.=vcovHC(mB2_h12_2,"HC2"))
mB2r_h13_1 <- coeftest(mB2_h13_1, vcov.=vcovHC(mB2_h13_1,"HC2"))
mB2r_h13_2 <- coeftest(mB2_h13_2, vcov.=vcovHC(mB2_h13_2,"HC2"))

# 仮説２の検証（３種類）
# 統制群vs実験群1＆実験群2vs実験群4 ＆ 実験群3vs実験群5
mB1_h20_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h20_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h21_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h21_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h22_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_self+1),ctl), data=dtmp)
mB1_h22_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_self+1),ctl), data=dtmp)

mB1r_h20_1 <- coeftest(mB1_h20_1, vcov.=vcovHC(mB1_h20_1,"HC2"))
mB1r_h20_2 <- coeftest(mB1_h20_2, vcov.=vcovHC(mB1_h20_2,"HC2"))
mB1r_h21_1 <- coeftest(mB1_h21_1, vcov.=vcovHC(mB1_h21_1,"HC2"))
mB1r_h21_2 <- coeftest(mB1_h21_2, vcov.=vcovHC(mB1_h21_2,"HC2"))
mB1r_h22_1 <- coeftest(mB1_h22_1, vcov.=vcovHC(mB1_h22_1,"HC2"))
mB1r_h22_2 <- coeftest(mB1_h22_2, vcov.=vcovHC(mB1_h22_2,"HC2"))

mB2_h20_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h20_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h21_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h21_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h22_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_self-1),ctl), data=dtmp)
mB2_h22_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_self-1),ctl), data=dtmp)

mB2r_h20_1 <- coeftest(mB2_h20_1, vcov.=vcovHC(mB2_h20_1,"HC2"))
mB2r_h20_2 <- coeftest(mB2_h20_2, vcov.=vcovHC(mB2_h20_2,"HC2"))
mB2r_h21_1 <- coeftest(mB2_h21_1, vcov.=vcovHC(mB2_h21_1,"HC2"))
mB2r_h21_2 <- coeftest(mB2_h21_2, vcov.=vcovHC(mB2_h21_2,"HC2"))
mB2r_h22_1 <- coeftest(mB2_h22_1, vcov.=vcovHC(mB2_h22_1,"HC2"))
mB2r_h22_2 <- coeftest(mB2_h22_2, vcov.=vcovHC(mB2_h22_2,"HC2"))

#'
#' ## 争点態度イデオロギー条件付け
#'

#'
#' ### 外交安全保障
#' 

mx_ctax1 <- lm(update(tax1_opi ~ as.factor(g_ctax_N)*ide_iss_1,ctl), data=dtmp)
coeftest(mx_ctax1, vcovHC(mx_ctax1, "HC2"))
mBa_1 <- mx_ctax1
mx_ctax2 <- lm(update(tax2_opi ~ as.factor(g_ctax_N)*ide_iss_1,ctl), data=dtmp)
coeftest(mx_ctax2, vcovHC(mx_ctax2, "HC2"))
mBa_2 <- mx_ctax2

# 仮説１の検証（４種類）
# 統制群vs実験群2＆実験群3vs実験群2＆実験群1vs実験群4＆実験群5vs実験群4
mBa1_h10_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h10_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h11_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h11_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h12_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h12_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h13_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h13_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_iss_1+1.35),ctl), data=dtmp)

mBa1r_h10_1 <- coeftest(mBa1_h10_1, vcov.=vcovHC(mBa1_h10_1,"HC2"))
mBa1r_h10_2 <- coeftest(mBa1_h10_2, vcov.=vcovHC(mBa1_h10_2,"HC2"))
mBa1r_h11_1 <- coeftest(mBa1_h11_1, vcov.=vcovHC(mBa1_h11_1,"HC2"))
mBa1r_h11_2 <- coeftest(mBa1_h11_2, vcov.=vcovHC(mBa1_h11_2,"HC2"))
mBa1r_h12_1 <- coeftest(mBa1_h12_1, vcov.=vcovHC(mBa1_h12_1,"HC2"))
mBa1r_h12_2 <- coeftest(mBa1_h12_2, vcov.=vcovHC(mBa1_h12_2,"HC2"))
mBa1r_h13_1 <- coeftest(mBa1_h13_1, vcov.=vcovHC(mBa1_h13_1,"HC2"))
mBa1r_h13_2 <- coeftest(mBa1_h13_2, vcov.=vcovHC(mBa1_h13_2,"HC2"))

mBa2_h10_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h10_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h11_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h11_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h12_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h12_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h13_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h13_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_iss_1-1.53),ctl), data=dtmp)

mBa2r_h10_1 <- coeftest(mBa2_h10_1, vcov.=vcovHC(mBa2_h10_1,"HC2"))
mBa2r_h10_2 <- coeftest(mBa2_h10_2, vcov.=vcovHC(mBa2_h10_2,"HC2"))
mBa2r_h11_1 <- coeftest(mBa2_h11_1, vcov.=vcovHC(mBa2_h11_1,"HC2"))
mBa2r_h11_2 <- coeftest(mBa2_h11_2, vcov.=vcovHC(mBa2_h11_2,"HC2"))
mBa2r_h12_1 <- coeftest(mBa2_h12_1, vcov.=vcovHC(mBa2_h12_1,"HC2"))
mBa2r_h12_2 <- coeftest(mBa2_h12_2, vcov.=vcovHC(mBa2_h12_2,"HC2"))
mBa2r_h13_1 <- coeftest(mBa2_h13_1, vcov.=vcovHC(mBa2_h13_1,"HC2"))
mBa2r_h13_2 <- coeftest(mBa2_h13_2, vcov.=vcovHC(mBa2_h13_2,"HC2"))

# 仮説２の検証（３種類）
# 統制群vs実験群1＆実験群2vs実験群4 ＆ 実験群3vs実験群5
mBa1_h20_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h20_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h21_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h21_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h22_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h22_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_iss_1+1.35),ctl), data=dtmp)

mBa1r_h20_1 <- coeftest(mBa1_h20_1, vcov.=vcovHC(mBa1_h20_1,"HC2"))
mBa1r_h20_2 <- coeftest(mBa1_h20_2, vcov.=vcovHC(mBa1_h20_2,"HC2"))
mBa1r_h21_1 <- coeftest(mBa1_h21_1, vcov.=vcovHC(mBa1_h21_1,"HC2"))
mBa1r_h21_2 <- coeftest(mBa1_h21_2, vcov.=vcovHC(mBa1_h21_2,"HC2"))
mBa1r_h22_1 <- coeftest(mBa1_h22_1, vcov.=vcovHC(mBa1_h22_1,"HC2"))
mBa1r_h22_2 <- coeftest(mBa1_h22_2, vcov.=vcovHC(mBa1_h22_2,"HC2"))

mBa2_h20_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h20_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h21_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h21_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h22_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h22_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_iss_1-1.53),ctl), data=dtmp)

mBa2r_h20_1 <- coeftest(mBa2_h20_1, vcov.=vcovHC(mBa2_h20_1,"HC2"))
mBa2r_h20_2 <- coeftest(mBa2_h20_2, vcov.=vcovHC(mBa2_h20_2,"HC2"))
mBa2r_h21_1 <- coeftest(mBa2_h21_1, vcov.=vcovHC(mBa2_h21_1,"HC2"))
mBa2r_h21_2 <- coeftest(mBa2_h21_2, vcov.=vcovHC(mBa2_h21_2,"HC2"))
mBa2r_h22_1 <- coeftest(mBa2_h22_1, vcov.=vcovHC(mBa2_h22_1,"HC2"))
mBa2r_h22_2 <- coeftest(mBa2_h22_2, vcov.=vcovHC(mBa2_h22_2,"HC2"))

#' 
#' ### 権利機会平等
#' 

mx_ctax1 <- lm(update(tax1_opi ~ as.factor(g_ctax_N)*ide_iss_2,ctl), data=dtmp)
coeftest(mx_ctax1, vcovHC(mx_ctax1, "HC2"))
mBb_1 <- mx_ctax1
mx_ctax2 <- lm(update(tax2_opi ~ as.factor(g_ctax_N)*ide_iss_2,ctl), data=dtmp)
coeftest(mx_ctax2, vcovHC(mx_ctax2, "HC2"))
mBb_2 <- mx_ctax2

# 仮説１の検証（４種類）
# 統制群vs実験群2＆実験群3vs実験群2＆実験群1vs実験群4＆実験群5vs実験群4
mBb1_h10_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h10_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h11_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h11_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h12_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h12_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h13_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h13_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_iss_2+1.50),ctl), data=dtmp)

mBb1r_h10_1 <- coeftest(mBb1_h10_1, vcov.=vcovHC(mBb1_h10_1,"HC2"))
mBb1r_h10_2 <- coeftest(mBb1_h10_2, vcov.=vcovHC(mBb1_h10_2,"HC2"))
mBb1r_h11_1 <- coeftest(mBb1_h11_1, vcov.=vcovHC(mBb1_h11_1,"HC2"))
mBb1r_h11_2 <- coeftest(mBb1_h11_2, vcov.=vcovHC(mBb1_h11_2,"HC2"))
mBb1r_h12_1 <- coeftest(mBb1_h12_1, vcov.=vcovHC(mBb1_h12_1,"HC2"))
mBb1r_h12_2 <- coeftest(mBb1_h12_2, vcov.=vcovHC(mBb1_h12_2,"HC2"))
mBb1r_h13_1 <- coeftest(mBb1_h13_1, vcov.=vcovHC(mBb1_h13_1,"HC2"))
mBb1r_h13_2 <- coeftest(mBb1_h13_2, vcov.=vcovHC(mBb1_h13_2,"HC2"))

mBb2_h10_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h10_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h11_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h11_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h12_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h12_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h13_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h13_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_iss_2-1.48),ctl), data=dtmp)

mBb2r_h10_1 <- coeftest(mBb2_h10_1, vcov.=vcovHC(mBb2_h10_1,"HC2"))
mBb2r_h10_2 <- coeftest(mBb2_h10_2, vcov.=vcovHC(mBb2_h10_2,"HC2"))
mBb2r_h11_1 <- coeftest(mBb2_h11_1, vcov.=vcovHC(mBb2_h11_1,"HC2"))
mBb2r_h11_2 <- coeftest(mBb2_h11_2, vcov.=vcovHC(mBb2_h11_2,"HC2"))
mBb2r_h12_1 <- coeftest(mBb2_h12_1, vcov.=vcovHC(mBb2_h12_1,"HC2"))
mBb2r_h12_2 <- coeftest(mBb2_h12_2, vcov.=vcovHC(mBb2_h12_2,"HC2"))
mBb2r_h13_1 <- coeftest(mBb2_h13_1, vcov.=vcovHC(mBb2_h13_1,"HC2"))
mBb2r_h13_2 <- coeftest(mBb2_h13_2, vcov.=vcovHC(mBb2_h13_2,"HC2"))

# 仮説２の検証（３種類）
# 統制群vs実験群1＆実験群2vs実験群4 ＆ 実験群3vs実験群5
mBb1_h20_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h20_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h21_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h21_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h22_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h22_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_iss_2+1.50),ctl), data=dtmp)

mBb1r_h20_1 <- coeftest(mBb1_h20_1, vcov.=vcovHC(mBb1_h20_1,"HC2"))
mBb1r_h20_2 <- coeftest(mBb1_h20_2, vcov.=vcovHC(mBb1_h20_2,"HC2"))
mBb1r_h21_1 <- coeftest(mBb1_h21_1, vcov.=vcovHC(mBb1_h21_1,"HC2"))
mBb1r_h21_2 <- coeftest(mBb1_h21_2, vcov.=vcovHC(mBb1_h21_2,"HC2"))
mBb1r_h22_1 <- coeftest(mBb1_h22_1, vcov.=vcovHC(mBb1_h22_1,"HC2"))
mBb1r_h22_2 <- coeftest(mBb1_h22_2, vcov.=vcovHC(mBb1_h22_2,"HC2"))

mBb2_h20_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h20_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h21_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h21_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h22_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h22_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_iss_2-1.48),ctl), data=dtmp)

mBb2r_h20_1 <- coeftest(mBb2_h20_1, vcov.=vcovHC(mBb2_h20_1,"HC2"))
mBb2r_h20_2 <- coeftest(mBb2_h20_2, vcov.=vcovHC(mBb2_h20_2,"HC2"))
mBb2r_h21_1 <- coeftest(mBb2_h21_1, vcov.=vcovHC(mBb2_h21_1,"HC2"))
mBb2r_h21_2 <- coeftest(mBb2_h21_2, vcov.=vcovHC(mBb2_h21_2,"HC2"))
mBb2r_h22_1 <- coeftest(mBb2_h22_1, vcov.=vcovHC(mBb2_h22_1,"HC2"))
mBb2r_h22_2 <- coeftest(mBb2_h22_2, vcov.=vcovHC(mBb2_h22_2,"HC2"))

#'
#' ## 政党支持イデオロギー条件付け
#'

mx_ctax1 <- lm(update(tax1_opi ~ as.factor(g_ctax_N)*ide_psup,ctl), data=dtmp)
coeftest(mx_ctax1, vcovHC(mx_ctax1, "HC2"))
mBc_1 <- mx_ctax1
mx_ctax2 <- lm(update(tax2_opi ~ as.factor(g_ctax_N)*ide_psup,ctl), data=dtmp)
coeftest(mx_ctax2, vcovHC(mx_ctax2, "HC2"))
mBc_2 <- mx_ctax2

# 仮説１の検証（４種類）
# 統制群vs実験群2＆実験群3vs実験群2＆実験群1vs実験群4＆実験群5vs実験群4
mBc1_h10_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h10_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h11_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h11_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h12_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h12_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h13_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h13_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_psup+1),ctl), data=dtmp)

mBc1r_h10_1 <- coeftest(mBc1_h10_1, vcov.=vcovHC(mBc1_h10_1,"HC2"))
mBc1r_h10_2 <- coeftest(mBc1_h10_2, vcov.=vcovHC(mBc1_h10_2,"HC2"))
mBc1r_h11_1 <- coeftest(mBc1_h11_1, vcov.=vcovHC(mBc1_h11_1,"HC2"))
mBc1r_h11_2 <- coeftest(mBc1_h11_2, vcov.=vcovHC(mBc1_h11_2,"HC2"))
mBc1r_h12_1 <- coeftest(mBc1_h12_1, vcov.=vcovHC(mBc1_h12_1,"HC2"))
mBc1r_h12_2 <- coeftest(mBc1_h12_2, vcov.=vcovHC(mBc1_h12_2,"HC2"))
mBc1r_h13_1 <- coeftest(mBc1_h13_1, vcov.=vcovHC(mBc1_h13_1,"HC2"))
mBc1r_h13_2 <- coeftest(mBc1_h13_2, vcov.=vcovHC(mBc1_h13_2,"HC2"))

mBc2_h10_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h10_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h11_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h11_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h12_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h12_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h13_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h13_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_psup-1),ctl), data=dtmp)

mBc2r_h10_1 <- coeftest(mBc2_h10_1, vcov.=vcovHC(mBc2_h10_1,"HC2"))
mBc2r_h10_2 <- coeftest(mBc2_h10_2, vcov.=vcovHC(mBc2_h10_2,"HC2"))
mBc2r_h11_1 <- coeftest(mBc2_h11_1, vcov.=vcovHC(mBc2_h11_1,"HC2"))
mBc2r_h11_2 <- coeftest(mBc2_h11_2, vcov.=vcovHC(mBc2_h11_2,"HC2"))
mBc2r_h12_1 <- coeftest(mBc2_h12_1, vcov.=vcovHC(mBc2_h12_1,"HC2"))
mBc2r_h12_2 <- coeftest(mBc2_h12_2, vcov.=vcovHC(mBc2_h12_2,"HC2"))
mBc2r_h13_1 <- coeftest(mBc2_h13_1, vcov.=vcovHC(mBc2_h13_1,"HC2"))
mBc2r_h13_2 <- coeftest(mBc2_h13_2, vcov.=vcovHC(mBc2_h13_2,"HC2"))

# 仮説２の検証（３種類）
# 統制群vs実験群1＆実験群2vs実験群4 ＆ 実験群3vs実験群5
mBc1_h20_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h20_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h21_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h21_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h22_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h22_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_psup+1),ctl), data=dtmp)

mBc1r_h20_1 <- coeftest(mBc1_h20_1, vcov.=vcovHC(mBc1_h20_1,"HC2"))
mBc1r_h20_2 <- coeftest(mBc1_h20_2, vcov.=vcovHC(mBc1_h20_2,"HC2"))
mBc1r_h21_1 <- coeftest(mBc1_h21_1, vcov.=vcovHC(mBc1_h21_1,"HC2"))
mBc1r_h21_2 <- coeftest(mBc1_h21_2, vcov.=vcovHC(mBc1_h21_2,"HC2"))
mBc1r_h22_1 <- coeftest(mBc1_h22_1, vcov.=vcovHC(mBc1_h22_1,"HC2"))
mBc1r_h22_2 <- coeftest(mBc1_h22_2, vcov.=vcovHC(mBc1_h22_2,"HC2"))

mBc2_h20_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h20_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h21_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h21_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h22_1 <- lm(update(tax1_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h22_2 <- lm(update(tax2_opi ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_psup-1),ctl), data=dtmp)

mBc2r_h20_1 <- coeftest(mBc2_h20_1, vcov.=vcovHC(mBc2_h20_1,"HC2"))
mBc2r_h20_2 <- coeftest(mBc2_h20_2, vcov.=vcovHC(mBc2_h20_2,"HC2"))
mBc2r_h21_1 <- coeftest(mBc2_h21_1, vcov.=vcovHC(mBc2_h21_1,"HC2"))
mBc2r_h21_2 <- coeftest(mBc2_h21_2, vcov.=vcovHC(mBc2_h21_2,"HC2"))
mBc2r_h22_1 <- coeftest(mBc2_h22_1, vcov.=vcovHC(mBc2_h22_1,"HC2"))
mBc2r_h22_2 <- coeftest(mBc2_h22_2, vcov.=vcovHC(mBc2_h22_2,"HC2"))

#'
#' ## 交差項による仮説検証
#'   

htest0 <- data.frame(int = rep(c("世帯収入", "自己申告イデオロギー",
                                 "外交安全保障イデオロギー",
                                 "権利機会平等イデオロギー",
                                 "政党支持イデオロギー"
), each=14),
dv = rep(c("生活必需品","その他すべて"), each=7),
h = rep(c("H1A/B","H1A/B",
          "H1A/B","H1A/B",
          "H2A/B","H2A/B",
          "H2A/B"),2*5),
cp = rep(c("2.普遍 - 0.統制",
           "2.普遍 - 3.選別",
           "4.逆進+普遍 - 1.逆進",
           "4.逆進+普遍 - 5.逆進+選別",
           "1.逆進 - 0.統制",
           "4.普遍+逆進 - 2.普遍",
           "5.選別+逆進 - 3.選別"),2*5),
rbind(mA1r_h10_1[17,],mA1r_h11_1[17,],mA1r_h12_1[17,],mA1r_h13_1[17,],
      mA1r_h20_1[17,],mA1r_h21_1[17,],mA1r_h22_1[17,],
      mA1r_h10_2[17,],mA1r_h11_2[17,],mA1r_h12_2[17,],mA1r_h13_2[17,],
      mA1r_h20_2[17,],mA1r_h21_2[17,],mA1r_h22_2[17,],
      mB1r_h10_1[17,],mB1r_h11_1[17,],mB1r_h12_1[17,],mB1r_h13_1[17,],
      mB1r_h20_1[17,],mB1r_h21_1[17,],mB1r_h22_1[17,],
      mB1r_h10_2[17,],mB1r_h11_2[17,],mB1r_h12_2[17,],mB1r_h13_2[17,],
      mB1r_h20_2[17,],mB1r_h21_2[17,],mB1r_h22_2[17,],
      mBa1r_h10_1[17,],mBa1r_h11_1[17,],mBa1r_h12_1[17,],mBa1r_h13_1[17,],
      mBa1r_h20_1[17,],mBa1r_h21_1[17,],mBa1r_h22_1[17,],
      mBa1r_h10_2[17,],mBa1r_h11_2[17,],mBa1r_h12_2[17,],mBa1r_h13_2[17,],
      mBa1r_h20_2[17,],mBa1r_h21_2[17,],mBa1r_h22_2[17,],
      mBb1r_h10_1[17,],mBb1r_h11_1[17,],mBb1r_h12_1[17,],mBb1r_h13_1[17,],
      mBb1r_h20_1[17,],mBb1r_h21_1[17,],mBb1r_h22_1[17,],
      mBb1r_h10_2[17,],mBb1r_h11_2[17,],mBb1r_h12_2[17,],mBb1r_h13_2[17,],
      mBb1r_h20_2[17,],mBb1r_h21_2[17,],mBb1r_h22_2[17,],
      mBc1r_h10_1[17,],mBc1r_h11_1[17,],mBc1r_h12_1[17,],mBc1r_h13_1[17,],
      mBc1r_h20_1[17,],mBc1r_h21_1[17,],mBc1r_h22_1[17,],
      mBc1r_h10_2[17,],mBc1r_h11_2[17,],mBc1r_h12_2[17,],mBc1r_h13_2[17,],
      mBc1r_h20_2[17,],mBc1r_h21_2[17,],mBc1r_h22_2[17,]))
htest0$int <- factor(htest0$int, levels=unique(htest0$int))
htest0$dv <- factor(htest0$dv, levels=unique(htest0$dv))
htest0$cp <- factor(htest0$cp, levels=rev(unique(htest0$cp)))
htest0$lo95 <- htest0$Estimate - qnorm(0.975)*htest0$Std..Error
htest0$up95 <- htest0$Estimate + qnorm(0.975)*htest0$Std..Error
htest0$lo90 <- htest0$Estimate - qnorm(0.95)*htest0$Std..Error
htest0$up90 <- htest0$Estimate + qnorm(0.95)*htest0$Std..Error

#' 
#' ### 実験情報刺激効果と収入・イデオロギーの交差項係数による仮説検証（その他のイデオロギー指標を含む）（図A13）
#'   

p <- ggplot(subset(htest0, int%in%c("世帯収入","自己申告イデオロギー")), 
            aes(x=cp, y=Estimate)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lo95,ymax=up95, color=int, alpha=int), 
                width=0.25, position = position_dodge(width=-0.7)) + 
  geom_errorbar(aes(ymin=lo90,ymax=up90, color=int, alpha=int), 
                width=0, size=1.5, position = position_dodge(width=-0.7)) + 
  geom_point(aes(color=int, shape=int, alpha=int), size=3, 
             position = position_dodge(width=-0.7)) +
  facet_grid(h~dv, scales = "free_y", space = "free_y") + 
  coord_flip() + 
  scale_color_manual(name="収入・イデオロギー", values=rep("black",2)) + 
  # scale_color_brewer(name="収入・イデオロギー", type="qual", palette=2) + 
  scale_shape_discrete(name="収入・イデオロギー") + 
  scale_alpha_manual(name = "収入・イデオロギー", values=c(1,rep(0.3,1))) + 
  labs(x=NULL, y="実験群比較変数と収入・イデオロギーの交差項係数\n（従属変数は理想消費税率）", 
       caption="分析の詳細は回帰表を参照。太線は90%信頼区間、細線は95%信頼区間を示している。",
       subtitle = "対象商品") + 
  theme_bw() + theme(legend.position="bottom",
                     plot.subtitle = element_text(hjust=0.5),
                     strip.placement = "outside")

#+ fig.width=8, fig.height=6, dev="png", dpi=300, echo=FALSE
p

#+ eval=FALSE
# ggsave("htest0_v1_originalscale_wokn.png", p, width=8, height=6)

#' 
#' ### 実験情報刺激効果と収入・イデオロギーの交差項係数による仮説検証（その他のイデオロギー指標を含む）（図A16）
#'   

#+ warning=FALSE, fig.width=9, fig.height=7, dev="png", dpi=300, echo=FALSE
p <- ggplot(htest0, 
            aes(x=cp, y=Estimate)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lo95,ymax=up95, color=int, alpha=int), 
                width=0.25, position = position_dodge(width=-0.7)) + 
  geom_errorbar(aes(ymin=lo90,ymax=up90, color=int, alpha=int), 
                width=0, size=1.5, position = position_dodge(width=-0.7)) + 
  geom_point(aes(color=int, shape=int, alpha=int), size=3, 
             position = position_dodge(width=-0.7)) +
  facet_grid(h~dv, scales = "free_y", space = "free_y") + 
  coord_flip() + 
  scale_color_manual(name="収入・イデオロギー", values=rep("black",5)) + 
  # scale_color_brewer(name="収入・イデオロギー", type="qual", palette=2) + 
  scale_shape_discrete(name="収入・イデオロギー") + 
  scale_alpha_manual(name = "収入・イデオロギー", values=c(1,rep(0.3,4))) + 
  labs(x=NULL, y="実験群比較変数と収入・イデオロギーの交差項係数\n（従属変数は理想消費税率）", 
       caption="分析の詳細は回帰表を参照。太線は90%信頼区間、細線は95%信頼区間を示している。",
       subtitle = "対象商品") + 
  theme_bw() + theme(legend.position="bottom",
                     plot.subtitle = element_text(hjust=0.5),
                     strip.placement = "outside") + 
  guides(color=guide_legend(nrow=2,byrow=TRUE),
         shape=guide_legend(nrow=2,byrow=TRUE),
         alpha=guide_legend(nrow=2,byrow=TRUE))
p

#+ warning=FALSE, eval=FALSE
ggsave("htest0_v2_originalscale_wokn.png", p, width=9, height=7)
#+

#'
#' ## 表のエクスポート
#' 
#' 
#' ### 直接効果
#' 

texreg(list(m0_1,m0_2), 
       override.se = list(coeftest(m0_1,vcovHC(m0_1,"HC2"))[,2],
                          coeftest(m0_2,vcovHC(m0_2,"HC2"))[,2]),
       override.pvalues = list(coeftest(m0_1,vcovHC(m0_1,"HC2"))[,4],
                               coeftest(m0_2,vcovHC(m0_2,"HC2"))[,4]),
       symbol = "+",
       single.row=TRUE, digits = 3, stars = c(0.001,0.01,0.05,0.1),
       custom.coef.map = vnmap,
       custom.model.names = c("1:生活必需品","2:その他すべて"),
       caption = "理想消費税率に実験情報刺激が与える効果（重回帰分析）",
       caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
       label="basetab_os_wokn", dcolumn = TRUE, booktabs = TRUE, use.packages = FALSE,
       custom.note = "%stars. 最小二乗法による重回帰分析、ロバスト標準誤差使用．")
# texreg(list(m0_1,m0_2), 
#        override.se = list(coeftest(m0_1,vcovHC(m0_1,"HC2"))[,2],
#                           coeftest(m0_2,vcovHC(m0_2,"HC2"))[,2]),
#        override.pvalues = list(coeftest(m0_1,vcovHC(m0_1,"HC2"))[,4],
#                                coeftest(m0_2,vcovHC(m0_2,"HC2"))[,4]),
#        # file = "../out_v4/basetab_originalscale_wokn.html", symbol = "&dagger;", 
#        file = "../out_v4/basetab_originalscale_wokn.tex", symbol = "\\dagger", 
#        single.row=TRUE, digits = 3, stars = c(0.001,0.01,0.05,0.1),
#        custom.coef.map = vnmap,
#        custom.model.names = c("1:生活必需品","2:その他すべて"),
#        caption = "理想消費税率に実験情報刺激が与える効果（重回帰分析）",
#        caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
#        label="basetab_os_wokn", dcolumn = TRUE, booktabs = TRUE, use.packages = FALSE,
#        custom.note = "%stars. 最小二乗法による重回帰分析、ロバスト標準誤差使用．")
# tmp <- readLines("../out_v4/basetab_originalscale_wokn.tex")
# tmp <- gsub("{dagger}","{\\dagger}", tmp, fixed=TRUE)
# writeLines(tmp, "../out_v4/basetab_originalscale_wokn.tex", useBytes = TRUE)

#' 
#' ### 収入・自己申告イデオロギーによる条件付け（表A4）
#' 

## Full Table
screenreg(list(mA_1,mA_2,mB_1,mB_2), 
          override.se = list(coeftest(mA_1,vcovHC(mA_1,"HC2"))[,2],
                             coeftest(mA_2,vcovHC(mA_2,"HC2"))[,2],
                             coeftest(mB_1,vcovHC(mB_1,"HC2"))[,2],
                             coeftest(mB_2,vcovHC(mB_2,"HC2"))[,2]),
          override.pvalues = list(coeftest(mA_1,vcovHC(mA_1,"HC2"))[,4],
                                  coeftest(mA_2,vcovHC(mA_2,"HC2"))[,4],
                                  coeftest(mB_1,vcovHC(mB_1,"HC2"))[,4],
                                  coeftest(mB_2,vcovHC(mB_2,"HC2"))[,4]),
          symbol = "+",
          single.row=FALSE, digits = 3, stars = c(0.001,0.01,0.05,0.1),
          custom.coef.map = vnmap,
          custom.model.names = c("1:生活必需品","2:その他すべて",
                                 "3:生活必需品","4:その他すべて"),
          custom.header = list("世帯収入" = 1:2, "自己申告イデオロギー" = 3:4), 
          caption = "理想消費税率に実験情報刺激が与える効果と収入・イデオロギー（重回帰分析）",
          caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
          label="maintab_full_os_wokn", dcolumn = TRUE, booktabs = TRUE, use.packages = FALSE,
          custom.note = "%stars. 最小二乗法による重回帰分析、ロバスト標準誤差使用．")
# texreg(list(mA_1,mA_2,mB_1,mB_2), 
#        override.se = list(coeftest(mA_1,vcovHC(mA_1,"HC2"))[,2],
#                           coeftest(mA_2,vcovHC(mA_2,"HC2"))[,2],
#                           coeftest(mB_1,vcovHC(mB_1,"HC2"))[,2],
#                           coeftest(mB_2,vcovHC(mB_2,"HC2"))[,2]),
#        override.pvalues = list(coeftest(mA_1,vcovHC(mA_1,"HC2"))[,4],
#                                coeftest(mA_2,vcovHC(mA_2,"HC2"))[,4],
#                                coeftest(mB_1,vcovHC(mB_1,"HC2"))[,4],
#                                coeftest(mB_2,vcovHC(mB_2,"HC2"))[,4]),
#        # file = "maintab_full_originalscale_wokn.html", symbol = "&dagger;", 
#        file = "maintab_full_originalscale_wokn.tex", symbol = "\\dagger", 
#        single.row=FALSE, digits = 3, stars = c(0.001,0.01,0.05,0.1),
#        custom.coef.map = vnmap,
#        custom.model.names = c("1:生活必需品","2:その他すべて",
#                               "3:生活必需品","4:その他すべて"),
#        custom.header = list("世帯収入" = 1:2, "自己申告イデオロギー" = 3:4), 
#        caption = "理想消費税率に実験情報刺激が与える効果と収入・イデオロギー（重回帰分析）",
#        caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
#        label="maintab_full_os_wokn", dcolumn = TRUE, booktabs = TRUE, use.packages = FALSE,
#        custom.note = "%stars. 最小二乗法による重回帰分析、ロバスト標準誤差使用．")
# tmp <- readLines("maintab_full_originalscale_wokn.tex")
# tmp <- gsub("{dagger}","{\\dagger}", tmp, fixed=TRUE)
# writeLines(tmp, "maintab_full_originalscale_wokn.tex", useBytes = TRUE)

#' 
#' ### その他のイデオロギーによる条件付け
#' 

## 他のイデオロギー（生活必需品）
screenreg(list(mBc_1,mBa_1,mBb_1), 
          override.se = list(coeftest(mBc_1,vcovHC(mBc_1,"HC2"))[,2],
                             coeftest(mBa_1,vcovHC(mBa_1,"HC2"))[,2],
                             coeftest(mBb_1,vcovHC(mBb_1,"HC2"))[,2]),
          override.pvalues = list(coeftest(mBc_1,vcovHC(mBc_1,"HC2"))[,4],
                                  coeftest(mBa_1,vcovHC(mBa_1,"HC2"))[,4],
                                  coeftest(mBb_1,vcovHC(mBb_1,"HC2"))[,4]),
          symbol="+",
          single.row=TRUE, digits = 3, stars = c(0.001,0.01,0.05,0.1),
          custom.coef.map = vnmap2,
          custom.model.names = c("政党支持","外交安全保障","権利機会平等"),
          caption = "生活必需品の理想消費税率に実験情報刺激が与える効果とイデオロギー（重回帰分析）",
          caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
          label="idetab1_os_wokn", dcolumn = TRUE, booktabs = TRUE, use.packages = FALSE,
          custom.note = "%stars. 最小二乗法による重回帰分析、ロバスト標準誤差使用．")
# texreg(list(mBc_1,mBa_1,mBb_1), 
#        override.se = list(coeftest(mBc_1,vcovHC(mBc_1,"HC2"))[,2],
#                           coeftest(mBa_1,vcovHC(mBa_1,"HC2"))[,2],
#                           coeftest(mBb_1,vcovHC(mBb_1,"HC2"))[,2]),
#        override.pvalues = list(coeftest(mBc_1,vcovHC(mBc_1,"HC2"))[,4],
#                                coeftest(mBa_1,vcovHC(mBa_1,"HC2"))[,4],
#                                coeftest(mBb_1,vcovHC(mBb_1,"HC2"))[,4]),
#        # file = "idetab1_originalscale_wokn.html", symbol = "&dagger;",
#        file = "idetab1_originalscale_wokn.tex", symbol = "\\dagger",
#        single.row=TRUE, digits = 3, stars = c(0.001,0.01,0.05,0.1),
#        custom.coef.map = vnmap2,
#        custom.model.names = c("政党支持","外交安全保障","権利機会平等"),
#        caption = "生活必需品の理想消費税率に実験情報刺激が与える効果とイデオロギー（重回帰分析）",
#        caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
#        label="idetab1_os_wokn", dcolumn = TRUE, booktabs = TRUE, use.packages = FALSE,
#        custom.note = "%stars. 最小二乗法による重回帰分析、ロバスト標準誤差使用．")
# tmp <- readLines("idetab1_originalscale_wokn.tex")
# tmp <- gsub("{dagger}","{\\dagger}", tmp, fixed=TRUE)
# writeLines(tmp, "idetab1_originalscale_wokn.tex", useBytes = TRUE)

## 他のイデオロギー（その他の商品）
screenreg(list(mBc_2,mBa_2,mBb_2), 
          override.se = list(coeftest(mBc_2,vcovHC(mBc_2,"HC2"))[,2],
                             coeftest(mBa_2,vcovHC(mBa_2,"HC2"))[,2],
                             coeftest(mBb_2,vcovHC(mBb_2,"HC2"))[,2]),
          override.pvalues = list(coeftest(mBc_2,vcovHC(mBc_2,"HC2"))[,4],
                                  coeftest(mBa_2,vcovHC(mBa_2,"HC2"))[,4],
                                  coeftest(mBb_2,vcovHC(mBb_2,"HC2"))[,4]),
          symbol = "+",
          single.row=TRUE, digits = 3, stars = c(0.001,0.01,0.05,0.1),
          custom.coef.map = vnmap2,
          custom.model.names = c("政党支持","外交安全保障","権利機会平等"),
          caption = "その他全ての商品の理想消費税率に実験情報刺激が与える効果とイデオロギー（重回帰分析）",
          caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
          label="idetab2_os_wokn", dcolumn = TRUE, booktabs = TRUE, use.packages = FALSE,
          custom.note = "%stars. 最小二乗法による重回帰分析、ロバスト標準誤差使用．")
# texreg(list(mBc_2,mBa_2,mBb_2), 
#        override.se = list(coeftest(mBc_2,vcovHC(mBc_2,"HC2"))[,2],
#                           coeftest(mBa_2,vcovHC(mBa_2,"HC2"))[,2],
#                           coeftest(mBb_2,vcovHC(mBb_2,"HC2"))[,2]),
#        override.pvalues = list(coeftest(mBc_2,vcovHC(mBc_2,"HC2"))[,4],
#                                coeftest(mBa_2,vcovHC(mBa_2,"HC2"))[,4],
#                                coeftest(mBb_2,vcovHC(mBb_2,"HC2"))[,4]),
#        # file = "idetab2_originalscale_wokn.html", symbol = "&dagger;",
#        file = "idetab2_originalscale_wokn.tex", symbol = "\\dagger",
#        single.row=TRUE, digits = 3, stars = c(0.001,0.01,0.05,0.1),
#        custom.coef.map = vnmap2,
#        custom.model.names = c("政党支持","外交安全保障","権利機会平等"),
#        caption = "その他全ての商品の理想消費税率に実験情報刺激が与える効果とイデオロギー（重回帰分析）",
#        caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
#        label="idetab2_os_wokn", dcolumn = TRUE, booktabs = TRUE, use.packages = FALSE,
#        custom.note = "%stars. 最小二乗法による重回帰分析、ロバスト標準誤差使用．")
# tmp <- readLines("idetab2_originalscale_wokn.tex")
# tmp <- gsub("{dagger}","{\\dagger}", tmp, fixed=TRUE)
# writeLines(tmp, "idetab2_originalscale_wokn.tex", useBytes = TRUE)

#'
#' # 実験群比較（平方根）
#' 
#' ## 準備

# 統制変数
ctl <- formula( ~ .+ knall + fem + age + lvlen + ownh + 
                  as.factor(edu3) + wk + mar + cld)

# 限界効果予測値計算用データ
preddata <- data.frame(g_ctax_N = seq(0,5,1))
preddata$knall = median(dtmp$knall,na.rm=TRUE)
preddata$fem = median(dtmp$fem, na.rm=TRUE)
preddata$age = median(dtmp$age, na.rm=TRUE)
preddata$lvlen = median(dtmp$lvlen, na.rm=TRUE)
preddata$ownh = median(dtmp$ownh, na.rm=TRUE)
preddata$edu3 = median(as.numeric(dtmp$edu3)-1, na.rm=TRUE)
preddata$wk = median(dtmp$wk, na.rm=TRUE)
preddata$mar = median(dtmp$mar, na.rm=TRUE)
preddata$cld = median(dtmp$cld, na.rm=TRUE)[1]

#'
#' ## 実験刺激の直接効果
#'

m_ctax1 <- lm(update(sqrt(tax1_opi) ~ as.factor(g_ctax_N),ctl), data=dtmp)
m0_1 <- m_ctax1
m_ctax2 <- lm(update(sqrt(tax2_opi) ~ as.factor(g_ctax_N),ctl), data=dtmp)
m0_2 <- m_ctax2

coeftest(m_ctax1, vcov.=vcovHC(m_ctax1,"HC2"))
coeftest(m_ctax2, vcov.=vcovHC(m_ctax2,"HC2"))

# 普遍性刺激効果の検証（４種類）
# 統制群vs実験群2＆実験群3vs実験群2＆実験群1vs実験群4＆実験群5vs実験群4
m0_h10_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5)),ctl), data=dtmp)
m0_h10_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5)),ctl), data=dtmp)
m0_h11_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5)),ctl), data=dtmp)
m0_h11_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5)),ctl), data=dtmp)
m0_h12_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5)),ctl), data=dtmp)
m0_h12_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5)),ctl), data=dtmp)
m0_h13_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3)),ctl), data=dtmp)
m0_h13_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3)),ctl), data=dtmp)

m0r_h10_1 <- coeftest(m0_h10_1, vcov.=vcovHC(m0_h10_1,"HC2"))
m0r_h10_2 <- coeftest(m0_h10_2, vcov.=vcovHC(m0_h10_2,"HC2"))
m0r_h11_1 <- coeftest(m0_h11_1, vcov.=vcovHC(m0_h11_1,"HC2"))
m0r_h11_2 <- coeftest(m0_h11_2, vcov.=vcovHC(m0_h11_2,"HC2"))
m0r_h12_1 <- coeftest(m0_h12_1, vcov.=vcovHC(m0_h12_1,"HC2"))
m0r_h12_2 <- coeftest(m0_h12_2, vcov.=vcovHC(m0_h12_2,"HC2"))
m0r_h13_1 <- coeftest(m0_h13_1, vcov.=vcovHC(m0_h13_1,"HC2"))
m0r_h13_2 <- coeftest(m0_h13_2, vcov.=vcovHC(m0_h13_2,"HC2"))

# 逆進性刺激効果の検証（３種類）
# 統制群vs実験群1＆実験群2vs実験群4 ＆ 実験群3vs実験群5
m0_h20_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5)),ctl), data=dtmp)
m0_h20_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5)),ctl), data=dtmp)
m0_h21_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5)),ctl), data=dtmp)
m0_h21_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5)),ctl), data=dtmp)
m0_h22_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4)),ctl), data=dtmp)
m0_h22_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4)),ctl), data=dtmp)

m0r_h20_1 <- coeftest(m0_h20_1, vcov.=vcovHC(m0_h20_1,"HC2"))
m0r_h20_2 <- coeftest(m0_h20_2, vcov.=vcovHC(m0_h20_2,"HC2"))
m0r_h21_1 <- coeftest(m0_h21_1, vcov.=vcovHC(m0_h21_1,"HC2"))
m0r_h21_2 <- coeftest(m0_h21_2, vcov.=vcovHC(m0_h21_2,"HC2"))
m0r_h22_1 <- coeftest(m0_h22_1, vcov.=vcovHC(m0_h22_1,"HC2"))
m0r_h22_2 <- coeftest(m0_h22_2, vcov.=vcovHC(m0_h22_2,"HC2"))

#' 
#' ### 仮説に関係する実験群比較に関連する直接効果（図A11）
#'

htest <- data.frame(dv = rep(c("生活必需品","その他すべて"), each=7),
                    h = rep(c("社会保障普遍性","社会保障普遍性",
                              "社会保障普遍性","社会保障普遍性",
                              "消費税逆進性","消費税逆進性",
                              "消費税逆進性"),2),
                    cp = rep(c("2.普遍 - 0.統制",
                               "2.普遍 - 3.選別",
                               "4.逆進+普遍 - 1.逆進",
                               "4.逆進+普遍 - 5.逆進+選別",
                               "1.逆進 - 0.統制",
                               "4.普遍+逆進 - 2.普遍",
                               "5.選別+逆進 - 3.選別"),2),
                    rbind(m0r_h10_1[2,],m0r_h11_1[2,],m0r_h12_1[2,],m0r_h13_1[2,],
                          m0r_h20_1[2,],m0r_h21_1[2,],m0r_h22_1[2,],
                          m0r_h10_2[2,],m0r_h11_2[2,],m0r_h12_2[2,],m0r_h13_2[2,],
                          m0r_h20_2[2,],m0r_h21_2[2,],m0r_h22_2[2,]))
htest$dv <- factor(htest$dv, levels=unique(htest$dv))
htest$cp <- factor(htest$cp, levels=rev(unique(htest$cp)))
htest$h <- factor(htest$h, levels=unique(htest$h))
htest$lo95 <- htest$Estimate - qnorm(0.975)*htest$Std..Error
htest$up95 <- htest$Estimate + qnorm(0.975)*htest$Std..Error
htest$lo90 <- htest$Estimate - qnorm(0.95)*htest$Std..Error
htest$up90 <- htest$Estimate + qnorm(0.95)*htest$Std..Error

p <- ggplot(htest, aes(x=cp, y=Estimate)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lo95,ymax=up95), width=0.25) + 
  geom_errorbar(aes(ymin=lo90,ymax=up90), width=0, size=1.5) + 
  geom_point(size=3) +
  facet_grid(h~dv, scale="free_y", switch="y") + 
  coord_flip() + 
  # scale_color_brewer(name="対象商品", type="qual", palette=2) + 
  labs(x=NULL, y="実験刺激効果（従属変数は理想消費税率の平方根）", 
       caption="分析の詳細は回帰表を参照。統制変数有。太線は90%信頼区間、細線は95%信頼区間を示している。", 
       subtitle = "対象商品") + 
  theme_bw() + 
  theme(plot.subtitle = element_text(hjust=0.5),
        strip.background.y = element_blank(),
        strip.text.y = element_blank(),
        strip.placement = "outside")

#+ warning=FALSE, fig.width=8, fig.height=5, dev="png", dpi=300, echo=FALSE
p

#+ warning=FALSE, eval=FALSE
# ggsave("htest_m0.png", p, width=8, height=5)

#'
#' ## 世帯収入条件付け
#'

mx_ctax1 <- lm(update(sqrt(tax1_opi) ~ as.factor(g_ctax_N)*inc,ctl), data=dtmp)
coeftest(mx_ctax1, vcovHC(mx_ctax1, "HC2"))
mA_1 <- mx_ctax1
mx_ctax2 <- lm(update(sqrt(tax2_opi) ~ as.factor(g_ctax_N)*inc,ctl), data=dtmp)
coeftest(mx_ctax2, vcovHC(mx_ctax2, "HC2"))
mA_2 <- mx_ctax2

# 仮説１の検証（４種類）
# 統制群vs実験群2＆実験群3vs実験群2＆実験群1vs実験群4＆実験群5vs実験群4
mA1_h10_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(inc-2),ctl), data=dtmp)
mA1_h10_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(inc-2),ctl), data=dtmp)
mA1_h11_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(inc-2),ctl), data=dtmp)
mA1_h11_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(inc-2),ctl), data=dtmp)
mA1_h12_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(inc-2),ctl), data=dtmp)
mA1_h12_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(inc-2),ctl), data=dtmp)
mA1_h13_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(inc-2),ctl), data=dtmp)
mA1_h13_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(inc-2),ctl), data=dtmp)

mA1r_h10_1 <- coeftest(mA1_h10_1, vcov.=vcovHC(mA1_h10_1,"HC2"))
mA1r_h10_2 <- coeftest(mA1_h10_2, vcov.=vcovHC(mA1_h10_2,"HC2"))
mA1r_h11_1 <- coeftest(mA1_h11_1, vcov.=vcovHC(mA1_h11_1,"HC2"))
mA1r_h11_2 <- coeftest(mA1_h11_2, vcov.=vcovHC(mA1_h11_2,"HC2"))
mA1r_h12_1 <- coeftest(mA1_h12_1, vcov.=vcovHC(mA1_h12_1,"HC2"))
mA1r_h12_2 <- coeftest(mA1_h12_2, vcov.=vcovHC(mA1_h12_2,"HC2"))
mA1r_h13_1 <- coeftest(mA1_h13_1, vcov.=vcovHC(mA1_h13_1,"HC2"))
mA1r_h13_2 <- coeftest(mA1_h13_2, vcov.=vcovHC(mA1_h13_2,"HC2"))

mA2_h10_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(inc-5),ctl), data=dtmp)
mA2_h10_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(inc-5),ctl), data=dtmp)
mA2_h11_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(inc-5),ctl), data=dtmp)
mA2_h11_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(inc-5),ctl), data=dtmp)
mA2_h12_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(inc-5),ctl), data=dtmp)
mA2_h12_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(inc-5),ctl), data=dtmp)
mA2_h13_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(inc-5),ctl), data=dtmp)
mA2_h13_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(inc-5),ctl), data=dtmp)

mA2r_h10_1 <- coeftest(mA2_h10_1, vcov.=vcovHC(mA2_h10_1,"HC2"))
mA2r_h10_2 <- coeftest(mA2_h10_2, vcov.=vcovHC(mA2_h10_2,"HC2"))
mA2r_h11_1 <- coeftest(mA2_h11_1, vcov.=vcovHC(mA2_h11_1,"HC2"))
mA2r_h11_2 <- coeftest(mA2_h11_2, vcov.=vcovHC(mA2_h11_2,"HC2"))
mA2r_h12_1 <- coeftest(mA2_h12_1, vcov.=vcovHC(mA2_h12_1,"HC2"))
mA2r_h12_2 <- coeftest(mA2_h12_2, vcov.=vcovHC(mA2_h12_2,"HC2"))
mA2r_h13_1 <- coeftest(mA2_h13_1, vcov.=vcovHC(mA2_h13_1,"HC2"))
mA2r_h13_2 <- coeftest(mA2_h13_2, vcov.=vcovHC(mA2_h13_2,"HC2"))

# 仮説２の検証（３種類）
# 統制群vs実験群1＆実験群2vs実験群4 ＆ 実験群3vs実験群5
mA1_h20_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(inc-2),ctl), data=dtmp)
mA1_h20_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(inc-2),ctl), data=dtmp)
mA1_h21_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(inc-2),ctl), data=dtmp)
mA1_h21_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(inc-2),ctl), data=dtmp)
mA1_h22_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(inc-2),ctl), data=dtmp)
mA1_h22_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(inc-2),ctl), data=dtmp)

mA1r_h20_1 <- coeftest(mA1_h20_1, vcov.=vcovHC(mA1_h20_1,"HC2"))
mA1r_h20_2 <- coeftest(mA1_h20_2, vcov.=vcovHC(mA1_h20_2,"HC2"))
mA1r_h21_1 <- coeftest(mA1_h21_1, vcov.=vcovHC(mA1_h21_1,"HC2"))
mA1r_h21_2 <- coeftest(mA1_h21_2, vcov.=vcovHC(mA1_h21_2,"HC2"))
mA1r_h22_1 <- coeftest(mA1_h22_1, vcov.=vcovHC(mA1_h22_1,"HC2"))
mA1r_h22_2 <- coeftest(mA1_h22_2, vcov.=vcovHC(mA1_h22_2,"HC2"))

mA2_h20_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(inc-5),ctl), data=dtmp)
mA2_h20_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(inc-5),ctl), data=dtmp)
mA2_h21_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(inc-5),ctl), data=dtmp)
mA2_h21_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(inc-5),ctl), data=dtmp)
mA2_h22_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(inc-5),ctl), data=dtmp)
mA2_h22_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(inc-5),ctl), data=dtmp)

mA2r_h20_1 <- coeftest(mA2_h20_1, vcov.=vcovHC(mA2_h20_1,"HC2"))
mA2r_h20_2 <- coeftest(mA2_h20_2, vcov.=vcovHC(mA2_h20_2,"HC2"))
mA2r_h21_1 <- coeftest(mA2_h21_1, vcov.=vcovHC(mA2_h21_1,"HC2"))
mA2r_h21_2 <- coeftest(mA2_h21_2, vcov.=vcovHC(mA2_h21_2,"HC2"))
mA2r_h22_1 <- coeftest(mA2_h22_1, vcov.=vcovHC(mA2_h22_1,"HC2"))
mA2r_h22_2 <- coeftest(mA2_h22_2, vcov.=vcovHC(mA2_h22_2,"HC2"))

#'
#' ### 世帯収入に条件付けされた実験情報刺激の効果係数を用いた仮説検証
#'

htest <- data.frame(int = rep(c("200〜400万円(10%=2)","800〜1000万円(90%=5)"), each=14), 
                    dv = rep(c("生活必需品","その他すべて"), each=7),
                    h = rep(c("H1A","H1A",
                              "H1A","H1A",
                              "H2A","H2A",
                              "H2A"),4),
                    cp = rep(c("2.普遍 - 0.統制",
                               "2.普遍 - 3.選別",
                               "4.逆進+普遍 - 1.逆進",
                               "4.逆進+普遍 - 5.逆進+選別",
                               "1.逆進 - 0.統制",
                               "4.普遍+逆進 - 2.普遍",
                               "5.選別+逆進 - 3.選別"),4),
                    rbind(mA1r_h10_1[2,],mA1r_h11_1[2,],mA1r_h12_1[2,],mA1r_h13_1[2,],
                          mA1r_h20_1[2,],mA1r_h21_1[2,],mA1r_h22_1[2,],
                          mA1r_h10_2[2,],mA1r_h11_2[2,],mA1r_h12_2[2,],mA1r_h13_2[2,],
                          mA1r_h20_2[2,],mA1r_h21_2[2,],mA1r_h22_2[2,],
                          mA2r_h10_1[2,],mA2r_h11_1[2,],mA2r_h12_1[2,],mA2r_h13_1[2,],
                          mA2r_h20_1[2,],mA2r_h21_1[2,],mA2r_h22_1[2,],
                          mA2r_h10_2[2,],mA2r_h11_2[2,],mA2r_h12_2[2,],mA2r_h13_2[2,],
                          mA2r_h20_2[2,],mA2r_h21_2[2,],mA2r_h22_2[2,]))
htest$dv <- factor(htest$dv, levels=unique(htest$dv))
htest$cp <- factor(htest$cp, levels=rev(unique(htest$cp)))
htest$lo95 <- htest$Estimate - qnorm(0.975)*htest$Std..Error
htest$up95 <- htest$Estimate + qnorm(0.975)*htest$Std..Error
htest$lo90 <- htest$Estimate - qnorm(0.95)*htest$Std..Error
htest$up90 <- htest$Estimate + qnorm(0.95)*htest$Std..Error

#+ warning=FALSE, fig.width=8, fig.height=5, dev="png", dpi=300, echo=FALSE
p <- ggplot(htest, aes(x=cp, y=Estimate)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lo95,ymax=up95, color=int), width=0.25, position = position_dodge(width=-0.5)) + 
  geom_errorbar(aes(ymin=lo90,ymax=up90, color=int), width=0, size=1.5, position = position_dodge(width=-0.5)) + 
  geom_point(aes(color=int, shape=int), size=3, position = position_dodge(width=-0.5)) +
  facet_grid(h~dv, scales = "free_y", space = "free_y") + 
  coord_flip() + 
  scale_color_brewer(name="世帯収入", type="qual", palette=2) + 
  scale_shape_discrete(name="世帯収入") + 
  labs(x=NULL, y="実験刺激効果（従属変数は理想消費税率の平方根）", 
       caption="分析の詳細は回帰表を参照。太線は90%信頼区間、細線は95%信頼区間を示している。",
       subtitle = "対象商品") + 
  theme_bw() + theme(legend.position="bottom",
                     plot.subtitle = element_text(hjust=0.5),
                     strip.placement = "outside")
p

#+ warning=FALSE, eval=FALSE
ggsave("htest_mA.png", p, width=8, height=6)
#+

#'
#' ### 世帯収入に条件付けされた実験情報刺激の限界効果を用いた仮説検証（図A10）
#'

### Data for Simulation
prdt1 <- data.frame(preddata[rep(1,14),], inc=rep(c(2,5),each=7), tax1_opi=0, tax2_opi=0)
prdt1$g_ctax_N <- c(0,3,1,5,0,2,3)
prdt2 <- data.frame(preddata[rep(1,14),], inc=rep(c(2,5),each=7), tax1_opi=0, tax2_opi=0)
prdt2$g_ctax_N <- c(2,2,4,4,1,4,5)

### Set random number seed
set.seed(500)
### Draw Coefficients 1000 times
require(MASS)
m <- mA_1
betadraw <- mvrnorm(1000, m$coefficients, vcovHC(m, "HC2")) 
### Prepare Storage of Output
hmargin_1<-matrix(NA,nrow=nrow(prdt1),ncol=7) 
colnames(hmargin_1)<-c("Mean","Median","SE","lCI95","lCI90","uCI90","uCI95")
### Simulation
for(i in 1:nrow(prdt1)){
  ## 1st Profile (treated)
  xfix1 <- model.matrix(m, data=prdt1)[i,] 
  ## 2nd Profile (control)
  xfix2 <- model.matrix(m, data=prdt2)[i,]
  ## Predicted Probabilities
  pred1 <- (betadraw%*%xfix1)^2 ## for 1st profile
  pred2 <- (betadraw%*%xfix2)^2 ## for 2nd profile
  ## Difference in Predicted Values
  predstore <- pred2 - pred1
  ## Store Stats
  hmargin_1[i,]<-c(mean(predstore),median(predstore),sd(predstore),
                   quantile(predstore, probs=c(0.025,0.05,0.95,0.975))) # 95% CI
  ## Print progress
  if (i==1) tmpper <- -1
  if (round(i/nrow(prdt1)*100)>tmpper) {
    tmpper <- round(i/nrow(prdt1)*100)
    cat(paste0(tmpper,"% ")) 
    Sys.sleep(0.01)
  }
}

### Set random number seed
set.seed(500)
### Draw Coefficients 1000 times
require(MASS)
m <- mA_2
betadraw <- mvrnorm(1000, m$coefficients, vcovHC(m, "HC2")) 
### Prepare Storage of Output
hmargin_2<-matrix(NA,nrow=nrow(prdt1),ncol=7) 
colnames(hmargin_2)<-c("Mean","Median","SE","lCI95","lCI90","uCI90","uCI95")
### Simulation
for(i in 1:nrow(prdt1)){
  ## 1st Profile (Female without Kid)
  xfix1 <- model.matrix(m, data=prdt1)[i,] 
  ## 2nd Profile (Female with Kid)
  xfix2 <- model.matrix(m, data=prdt2)[i,]
  ## Predicted Probabilities
  pred1 <- (betadraw%*%xfix1)^2 ## for 1st profile
  pred2 <- (betadraw%*%xfix2)^2 ## for 2nd profile
  ## Difference in Predicted Values
  predstore <- pred2 - pred1
  ## Store Stats
  hmargin_2[i,]<-c(mean(predstore),median(predstore),sd(predstore),
                   quantile(predstore, probs=c(0.025,0.05,0.95,0.975))) # 95% CI
  ## Print progress
  if (i==1) tmpper <- -1
  if (round(i/nrow(prdt1)*100)>tmpper) {
    tmpper <- round(i/nrow(prdt1)*100)
    cat(paste0(tmpper,"% ")) 
    Sys.sleep(0.01)
  }
}

### Combine All the Information
hmargin <- cbind(rbind(hmargin_1, hmargin_2), htest[c(1:7,15:21,8:14,22:28),c(1:4)])

p <- ggplot(hmargin, aes(x=cp, y=Mean)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lCI95,ymax=uCI95, color=int), width=0.25, position = position_dodge(width=-0.5)) + 
  geom_errorbar(aes(ymin=lCI90,ymax=uCI90, color=int), width=0, size=1.5, position = position_dodge(width=-0.5)) + 
  geom_point(aes(color=int, shape=int), size=3, position = position_dodge(width=-0.5)) +
  facet_grid(h~dv, scales = "free_y", space = "free_y") + 
  coord_flip() + 
  scale_color_brewer(name="世帯収入", type="qual", palette=2) + 
  scale_shape_discrete(name="世帯収入") + 
  labs(x=NULL, y="実験刺激の限界効果（理想消費税率の変化量の平均値）", 
       caption="分析の詳細は回帰表を参照。太線は90%、細線は95%信頼区間を示している。統制変数を中央値で固定し、モンテカルロ・シミュレーションで推定。",
       subtitle = "対象商品") + 
  theme_bw() + theme(legend.position="bottom",
                     plot.subtitle = element_text(hjust=0.5),
                     strip.placement = "outside")

#+ warning=FALSE, fig.width=9, fig.height=6, dev="png", dpi=300, echo=FALSE
p

#+ warning=FALSE, eval=FALSE
# ggsave("hmargin_mA.png", p, width=9, height=6)

#'
#' ## 自己申告イデオロギー条件付け
#'

mx_ctax1 <- lm(update(sqrt(tax1_opi) ~ as.factor(g_ctax_N)*ide_self,ctl), data=dtmp)
coeftest(mx_ctax1, vcovHC(mx_ctax1, "HC2"))
mB_1 <- mx_ctax1
mx_ctax2 <- lm(update(sqrt(tax2_opi) ~ as.factor(g_ctax_N)*ide_self,ctl), data=dtmp)
coeftest(mx_ctax2, vcovHC(mx_ctax2, "HC2"))
mB_2 <- mx_ctax2

# 仮説１の検証（４種類）
# 統制群vs実験群2＆実験群3vs実験群2＆実験群1vs実験群4＆実験群5vs実験群4
mB1_h10_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h10_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h11_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h11_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h12_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h12_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h13_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_self+1),ctl), data=dtmp)
mB1_h13_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_self+1),ctl), data=dtmp)

mB1r_h10_1 <- coeftest(mB1_h10_1, vcov.=vcovHC(mB1_h10_1,"HC2"))
mB1r_h10_2 <- coeftest(mB1_h10_2, vcov.=vcovHC(mB1_h10_2,"HC2"))
mB1r_h11_1 <- coeftest(mB1_h11_1, vcov.=vcovHC(mB1_h11_1,"HC2"))
mB1r_h11_2 <- coeftest(mB1_h11_2, vcov.=vcovHC(mB1_h11_2,"HC2"))
mB1r_h12_1 <- coeftest(mB1_h12_1, vcov.=vcovHC(mB1_h12_1,"HC2"))
mB1r_h12_2 <- coeftest(mB1_h12_2, vcov.=vcovHC(mB1_h12_2,"HC2"))
mB1r_h13_1 <- coeftest(mB1_h13_1, vcov.=vcovHC(mB1_h13_1,"HC2"))
mB1r_h13_2 <- coeftest(mB1_h13_2, vcov.=vcovHC(mB1_h13_2,"HC2"))

mB2_h10_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h10_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h11_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h11_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h12_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h12_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h13_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_self-1),ctl), data=dtmp)
mB2_h13_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_self-1),ctl), data=dtmp)

mB2r_h10_1 <- coeftest(mB2_h10_1, vcov.=vcovHC(mB2_h10_1,"HC2"))
mB2r_h10_2 <- coeftest(mB2_h10_2, vcov.=vcovHC(mB2_h10_2,"HC2"))
mB2r_h11_1 <- coeftest(mB2_h11_1, vcov.=vcovHC(mB2_h11_1,"HC2"))
mB2r_h11_2 <- coeftest(mB2_h11_2, vcov.=vcovHC(mB2_h11_2,"HC2"))
mB2r_h12_1 <- coeftest(mB2_h12_1, vcov.=vcovHC(mB2_h12_1,"HC2"))
mB2r_h12_2 <- coeftest(mB2_h12_2, vcov.=vcovHC(mB2_h12_2,"HC2"))
mB2r_h13_1 <- coeftest(mB2_h13_1, vcov.=vcovHC(mB2_h13_1,"HC2"))
mB2r_h13_2 <- coeftest(mB2_h13_2, vcov.=vcovHC(mB2_h13_2,"HC2"))

# 仮説２の検証（３種類）
# 統制群vs実験群1＆実験群2vs実験群4 ＆ 実験群3vs実験群5
mB1_h20_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h20_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h21_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h21_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_self+1),ctl), data=dtmp)
mB1_h22_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_self+1),ctl), data=dtmp)
mB1_h22_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_self+1),ctl), data=dtmp)

mB1r_h20_1 <- coeftest(mB1_h20_1, vcov.=vcovHC(mB1_h20_1,"HC2"))
mB1r_h20_2 <- coeftest(mB1_h20_2, vcov.=vcovHC(mB1_h20_2,"HC2"))
mB1r_h21_1 <- coeftest(mB1_h21_1, vcov.=vcovHC(mB1_h21_1,"HC2"))
mB1r_h21_2 <- coeftest(mB1_h21_2, vcov.=vcovHC(mB1_h21_2,"HC2"))
mB1r_h22_1 <- coeftest(mB1_h22_1, vcov.=vcovHC(mB1_h22_1,"HC2"))
mB1r_h22_2 <- coeftest(mB1_h22_2, vcov.=vcovHC(mB1_h22_2,"HC2"))

mB2_h20_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h20_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h21_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h21_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_self-1),ctl), data=dtmp)
mB2_h22_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_self-1),ctl), data=dtmp)
mB2_h22_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_self-1),ctl), data=dtmp)

mB2r_h20_1 <- coeftest(mB2_h20_1, vcov.=vcovHC(mB2_h20_1,"HC2"))
mB2r_h20_2 <- coeftest(mB2_h20_2, vcov.=vcovHC(mB2_h20_2,"HC2"))
mB2r_h21_1 <- coeftest(mB2_h21_1, vcov.=vcovHC(mB2_h21_1,"HC2"))
mB2r_h21_2 <- coeftest(mB2_h21_2, vcov.=vcovHC(mB2_h21_2,"HC2"))
mB2r_h22_1 <- coeftest(mB2_h22_1, vcov.=vcovHC(mB2_h22_1,"HC2"))
mB2r_h22_2 <- coeftest(mB2_h22_2, vcov.=vcovHC(mB2_h22_2,"HC2"))

#'
#' ## 争点態度イデオロギー条件付け
#'

#'
#' ### 外交安全保障
#' 
mx_ctax1 <- lm(update(sqrt(tax1_opi) ~ as.factor(g_ctax_N)*ide_iss_1,ctl), data=dtmp)
coeftest(mx_ctax1, vcovHC(mx_ctax1, "HC2"))
mBa_1 <- mx_ctax1
mx_ctax2 <- lm(update(sqrt(tax2_opi) ~ as.factor(g_ctax_N)*ide_iss_1,ctl), data=dtmp)
coeftest(mx_ctax2, vcovHC(mx_ctax2, "HC2"))
mBa_2 <- mx_ctax2

# 仮説１の検証（４種類）
# 統制群vs実験群2＆実験群3vs実験群2＆実験群1vs実験群4＆実験群5vs実験群4
mBa1_h10_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h10_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h11_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h11_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h12_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h12_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h13_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h13_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_iss_1+1.35),ctl), data=dtmp)

mBa1r_h10_1 <- coeftest(mBa1_h10_1, vcov.=vcovHC(mBa1_h10_1,"HC2"))
mBa1r_h10_2 <- coeftest(mBa1_h10_2, vcov.=vcovHC(mBa1_h10_2,"HC2"))
mBa1r_h11_1 <- coeftest(mBa1_h11_1, vcov.=vcovHC(mBa1_h11_1,"HC2"))
mBa1r_h11_2 <- coeftest(mBa1_h11_2, vcov.=vcovHC(mBa1_h11_2,"HC2"))
mBa1r_h12_1 <- coeftest(mBa1_h12_1, vcov.=vcovHC(mBa1_h12_1,"HC2"))
mBa1r_h12_2 <- coeftest(mBa1_h12_2, vcov.=vcovHC(mBa1_h12_2,"HC2"))
mBa1r_h13_1 <- coeftest(mBa1_h13_1, vcov.=vcovHC(mBa1_h13_1,"HC2"))
mBa1r_h13_2 <- coeftest(mBa1_h13_2, vcov.=vcovHC(mBa1_h13_2,"HC2"))

mBa2_h10_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h10_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h11_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h11_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h12_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h12_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h13_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h13_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_iss_1-1.53),ctl), data=dtmp)

mBa2r_h10_1 <- coeftest(mBa2_h10_1, vcov.=vcovHC(mBa2_h10_1,"HC2"))
mBa2r_h10_2 <- coeftest(mBa2_h10_2, vcov.=vcovHC(mBa2_h10_2,"HC2"))
mBa2r_h11_1 <- coeftest(mBa2_h11_1, vcov.=vcovHC(mBa2_h11_1,"HC2"))
mBa2r_h11_2 <- coeftest(mBa2_h11_2, vcov.=vcovHC(mBa2_h11_2,"HC2"))
mBa2r_h12_1 <- coeftest(mBa2_h12_1, vcov.=vcovHC(mBa2_h12_1,"HC2"))
mBa2r_h12_2 <- coeftest(mBa2_h12_2, vcov.=vcovHC(mBa2_h12_2,"HC2"))
mBa2r_h13_1 <- coeftest(mBa2_h13_1, vcov.=vcovHC(mBa2_h13_1,"HC2"))
mBa2r_h13_2 <- coeftest(mBa2_h13_2, vcov.=vcovHC(mBa2_h13_2,"HC2"))

# 仮説２の検証（３種類）
# 統制群vs実験群1＆実験群2vs実験群4 ＆ 実験群3vs実験群5
mBa1_h20_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h20_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h21_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h21_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h22_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_iss_1+1.35),ctl), data=dtmp)
mBa1_h22_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_iss_1+1.35),ctl), data=dtmp)

mBa1r_h20_1 <- coeftest(mBa1_h20_1, vcov.=vcovHC(mBa1_h20_1,"HC2"))
mBa1r_h20_2 <- coeftest(mBa1_h20_2, vcov.=vcovHC(mBa1_h20_2,"HC2"))
mBa1r_h21_1 <- coeftest(mBa1_h21_1, vcov.=vcovHC(mBa1_h21_1,"HC2"))
mBa1r_h21_2 <- coeftest(mBa1_h21_2, vcov.=vcovHC(mBa1_h21_2,"HC2"))
mBa1r_h22_1 <- coeftest(mBa1_h22_1, vcov.=vcovHC(mBa1_h22_1,"HC2"))
mBa1r_h22_2 <- coeftest(mBa1_h22_2, vcov.=vcovHC(mBa1_h22_2,"HC2"))

mBa2_h20_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h20_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h21_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h21_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h22_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_iss_1-1.53),ctl), data=dtmp)
mBa2_h22_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_iss_1-1.53),ctl), data=dtmp)

mBa2r_h20_1 <- coeftest(mBa2_h20_1, vcov.=vcovHC(mBa2_h20_1,"HC2"))
mBa2r_h20_2 <- coeftest(mBa2_h20_2, vcov.=vcovHC(mBa2_h20_2,"HC2"))
mBa2r_h21_1 <- coeftest(mBa2_h21_1, vcov.=vcovHC(mBa2_h21_1,"HC2"))
mBa2r_h21_2 <- coeftest(mBa2_h21_2, vcov.=vcovHC(mBa2_h21_2,"HC2"))
mBa2r_h22_1 <- coeftest(mBa2_h22_1, vcov.=vcovHC(mBa2_h22_1,"HC2"))
mBa2r_h22_2 <- coeftest(mBa2_h22_2, vcov.=vcovHC(mBa2_h22_2,"HC2"))

#' 
#' ### 権利機会平等
#' 

mx_ctax1 <- lm(update(sqrt(tax1_opi) ~ as.factor(g_ctax_N)*ide_iss_2,ctl), data=dtmp)
coeftest(mx_ctax1, vcovHC(mx_ctax1, "HC2"))
mBb_1 <- mx_ctax1
mx_ctax2 <- lm(update(sqrt(tax2_opi) ~ as.factor(g_ctax_N)*ide_iss_2,ctl), data=dtmp)
coeftest(mx_ctax2, vcovHC(mx_ctax2, "HC2"))
mBb_2 <- mx_ctax2

# 仮説１の検証（４種類）
# 統制群vs実験群2＆実験群3vs実験群2＆実験群1vs実験群4＆実験群5vs実験群4
mBb1_h10_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h10_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h11_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h11_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h12_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h12_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h13_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h13_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_iss_2+1.50),ctl), data=dtmp)

mBb1r_h10_1 <- coeftest(mBb1_h10_1, vcov.=vcovHC(mBb1_h10_1,"HC2"))
mBb1r_h10_2 <- coeftest(mBb1_h10_2, vcov.=vcovHC(mBb1_h10_2,"HC2"))
mBb1r_h11_1 <- coeftest(mBb1_h11_1, vcov.=vcovHC(mBb1_h11_1,"HC2"))
mBb1r_h11_2 <- coeftest(mBb1_h11_2, vcov.=vcovHC(mBb1_h11_2,"HC2"))
mBb1r_h12_1 <- coeftest(mBb1_h12_1, vcov.=vcovHC(mBb1_h12_1,"HC2"))
mBb1r_h12_2 <- coeftest(mBb1_h12_2, vcov.=vcovHC(mBb1_h12_2,"HC2"))
mBb1r_h13_1 <- coeftest(mBb1_h13_1, vcov.=vcovHC(mBb1_h13_1,"HC2"))
mBb1r_h13_2 <- coeftest(mBb1_h13_2, vcov.=vcovHC(mBb1_h13_2,"HC2"))

mBb2_h10_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h10_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h11_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h11_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h12_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h12_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h13_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h13_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_iss_2-1.48),ctl), data=dtmp)

mBb2r_h10_1 <- coeftest(mBb2_h10_1, vcov.=vcovHC(mBb2_h10_1,"HC2"))
mBb2r_h10_2 <- coeftest(mBb2_h10_2, vcov.=vcovHC(mBb2_h10_2,"HC2"))
mBb2r_h11_1 <- coeftest(mBb2_h11_1, vcov.=vcovHC(mBb2_h11_1,"HC2"))
mBb2r_h11_2 <- coeftest(mBb2_h11_2, vcov.=vcovHC(mBb2_h11_2,"HC2"))
mBb2r_h12_1 <- coeftest(mBb2_h12_1, vcov.=vcovHC(mBb2_h12_1,"HC2"))
mBb2r_h12_2 <- coeftest(mBb2_h12_2, vcov.=vcovHC(mBb2_h12_2,"HC2"))
mBb2r_h13_1 <- coeftest(mBb2_h13_1, vcov.=vcovHC(mBb2_h13_1,"HC2"))
mBb2r_h13_2 <- coeftest(mBb2_h13_2, vcov.=vcovHC(mBb2_h13_2,"HC2"))

# 仮説２の検証（３種類）
# 統制群vs実験群1＆実験群2vs実験群4 ＆ 実験群3vs実験群5
mBb1_h20_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h20_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h21_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h21_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h22_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_iss_2+1.50),ctl), data=dtmp)
mBb1_h22_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_iss_2+1.50),ctl), data=dtmp)

mBb1r_h20_1 <- coeftest(mBb1_h20_1, vcov.=vcovHC(mBb1_h20_1,"HC2"))
mBb1r_h20_2 <- coeftest(mBb1_h20_2, vcov.=vcovHC(mBb1_h20_2,"HC2"))
mBb1r_h21_1 <- coeftest(mBb1_h21_1, vcov.=vcovHC(mBb1_h21_1,"HC2"))
mBb1r_h21_2 <- coeftest(mBb1_h21_2, vcov.=vcovHC(mBb1_h21_2,"HC2"))
mBb1r_h22_1 <- coeftest(mBb1_h22_1, vcov.=vcovHC(mBb1_h22_1,"HC2"))
mBb1r_h22_2 <- coeftest(mBb1_h22_2, vcov.=vcovHC(mBb1_h22_2,"HC2"))

mBb2_h20_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h20_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h21_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h21_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h22_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_iss_2-1.48),ctl), data=dtmp)
mBb2_h22_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_iss_2-1.48),ctl), data=dtmp)

mBb2r_h20_1 <- coeftest(mBb2_h20_1, vcov.=vcovHC(mBb2_h20_1,"HC2"))
mBb2r_h20_2 <- coeftest(mBb2_h20_2, vcov.=vcovHC(mBb2_h20_2,"HC2"))
mBb2r_h21_1 <- coeftest(mBb2_h21_1, vcov.=vcovHC(mBb2_h21_1,"HC2"))
mBb2r_h21_2 <- coeftest(mBb2_h21_2, vcov.=vcovHC(mBb2_h21_2,"HC2"))
mBb2r_h22_1 <- coeftest(mBb2_h22_1, vcov.=vcovHC(mBb2_h22_1,"HC2"))
mBb2r_h22_2 <- coeftest(mBb2_h22_2, vcov.=vcovHC(mBb2_h22_2,"HC2"))

#'
#' ## 政党支持イデオロギー条件付け
#'

mx_ctax1 <- lm(update(sqrt(tax1_opi) ~ as.factor(g_ctax_N)*ide_psup,ctl), data=dtmp)
coeftest(mx_ctax1, vcovHC(mx_ctax1, "HC2"))
mBc_1 <- mx_ctax1
mx_ctax2 <- lm(update(sqrt(tax2_opi) ~ as.factor(g_ctax_N)*ide_psup,ctl), data=dtmp)
coeftest(mx_ctax2, vcovHC(mx_ctax2, "HC2"))
mBc_2 <- mx_ctax2

# 仮説１の検証（４種類）
# 統制群vs実験群2＆実験群3vs実験群2＆実験群1vs実験群4＆実験群5vs実験群4
mBc1_h10_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h10_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h11_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h11_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h12_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h12_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h13_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h13_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_psup+1),ctl), data=dtmp)

mBc1r_h10_1 <- coeftest(mBc1_h10_1, vcov.=vcovHC(mBc1_h10_1,"HC2"))
mBc1r_h10_2 <- coeftest(mBc1_h10_2, vcov.=vcovHC(mBc1_h10_2,"HC2"))
mBc1r_h11_1 <- coeftest(mBc1_h11_1, vcov.=vcovHC(mBc1_h11_1,"HC2"))
mBc1r_h11_2 <- coeftest(mBc1_h11_2, vcov.=vcovHC(mBc1_h11_2,"HC2"))
mBc1r_h12_1 <- coeftest(mBc1_h12_1, vcov.=vcovHC(mBc1_h12_1,"HC2"))
mBc1r_h12_2 <- coeftest(mBc1_h12_2, vcov.=vcovHC(mBc1_h12_2,"HC2"))
mBc1r_h13_1 <- coeftest(mBc1_h13_1, vcov.=vcovHC(mBc1_h13_1,"HC2"))
mBc1r_h13_2 <- coeftest(mBc1_h13_2, vcov.=vcovHC(mBc1_h13_2,"HC2"))

mBc2_h10_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h10_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(0,2,1,3,4,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h11_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h11_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(3,2,0,1,4,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h12_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h12_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(1,4,0,2,3,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h13_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h13_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(5,4,0,1,2,3))*I(ide_psup-1),ctl), data=dtmp)

mBc2r_h10_1 <- coeftest(mBc2_h10_1, vcov.=vcovHC(mBc2_h10_1,"HC2"))
mBc2r_h10_2 <- coeftest(mBc2_h10_2, vcov.=vcovHC(mBc2_h10_2,"HC2"))
mBc2r_h11_1 <- coeftest(mBc2_h11_1, vcov.=vcovHC(mBc2_h11_1,"HC2"))
mBc2r_h11_2 <- coeftest(mBc2_h11_2, vcov.=vcovHC(mBc2_h11_2,"HC2"))
mBc2r_h12_1 <- coeftest(mBc2_h12_1, vcov.=vcovHC(mBc2_h12_1,"HC2"))
mBc2r_h12_2 <- coeftest(mBc2_h12_2, vcov.=vcovHC(mBc2_h12_2,"HC2"))
mBc2r_h13_1 <- coeftest(mBc2_h13_1, vcov.=vcovHC(mBc2_h13_1,"HC2"))
mBc2r_h13_2 <- coeftest(mBc2_h13_2, vcov.=vcovHC(mBc2_h13_2,"HC2"))

# 仮説２の検証（３種類）
# 統制群vs実験群1＆実験群2vs実験群4 ＆ 実験群3vs実験群5
mBc1_h20_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h20_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h21_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h21_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h22_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_psup+1),ctl), data=dtmp)
mBc1_h22_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_psup+1),ctl), data=dtmp)

mBc1r_h20_1 <- coeftest(mBc1_h20_1, vcov.=vcovHC(mBc1_h20_1,"HC2"))
mBc1r_h20_2 <- coeftest(mBc1_h20_2, vcov.=vcovHC(mBc1_h20_2,"HC2"))
mBc1r_h21_1 <- coeftest(mBc1_h21_1, vcov.=vcovHC(mBc1_h21_1,"HC2"))
mBc1r_h21_2 <- coeftest(mBc1_h21_2, vcov.=vcovHC(mBc1_h21_2,"HC2"))
mBc1r_h22_1 <- coeftest(mBc1_h22_1, vcov.=vcovHC(mBc1_h22_1,"HC2"))
mBc1r_h22_2 <- coeftest(mBc1_h22_2, vcov.=vcovHC(mBc1_h22_2,"HC2"))

mBc2_h20_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h20_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(0,1,2,3,4,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h21_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h21_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(2,4,0,1,3,5))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h22_1 <- lm(update(sqrt(tax1_opi) ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_psup-1),ctl), data=dtmp)
mBc2_h22_2 <- lm(update(sqrt(tax2_opi) ~ factor(g_ctax_N, levels=c(3,5,0,1,2,4))*I(ide_psup-1),ctl), data=dtmp)

mBc2r_h20_1 <- coeftest(mBc2_h20_1, vcov.=vcovHC(mBc2_h20_1,"HC2"))
mBc2r_h20_2 <- coeftest(mBc2_h20_2, vcov.=vcovHC(mBc2_h20_2,"HC2"))
mBc2r_h21_1 <- coeftest(mBc2_h21_1, vcov.=vcovHC(mBc2_h21_1,"HC2"))
mBc2r_h21_2 <- coeftest(mBc2_h21_2, vcov.=vcovHC(mBc2_h21_2,"HC2"))
mBc2r_h22_1 <- coeftest(mBc2_h22_1, vcov.=vcovHC(mBc2_h22_1,"HC2"))
mBc2r_h22_2 <- coeftest(mBc2_h22_2, vcov.=vcovHC(mBc2_h22_2,"HC2"))

#'
#' ## 交差項による仮説検証
#'   

htest0 <- data.frame(int = rep(c("世帯収入", "自己申告イデオロギー",
                                 "外交安全保障イデオロギー",
                                 "権利機会平等イデオロギー",
                                 "政党支持イデオロギー"
), each=14),
dv = rep(c("生活必需品","その他すべて"), each=7),
h = rep(c("H1A/B","H1A/B",
          "H1A/B","H1A/B",
          "H2A/B","H2A/B",
          "H2A/B"),2*5),
cp = rep(c("2.普遍 - 0.統制",
           "2.普遍 - 3.選別",
           "4.逆進+普遍 - 1.逆進",
           "4.逆進+普遍 - 5.逆進+選別",
           "1.逆進 - 0.統制",
           "4.普遍+逆進 - 2.普遍",
           "5.選別+逆進 - 3.選別"),2*5),
rbind(mA1r_h10_1[18,],mA1r_h11_1[18,],mA1r_h12_1[18,],mA1r_h13_1[18,],
      mA1r_h20_1[18,],mA1r_h21_1[18,],mA1r_h22_1[18,],
      mA1r_h10_2[18,],mA1r_h11_2[18,],mA1r_h12_2[18,],mA1r_h13_2[18,],
      mA1r_h20_2[18,],mA1r_h21_2[18,],mA1r_h22_2[18,],
      mB1r_h10_1[18,],mB1r_h11_1[18,],mB1r_h12_1[18,],mB1r_h13_1[18,],
      mB1r_h20_1[18,],mB1r_h21_1[18,],mB1r_h22_1[18,],
      mB1r_h10_2[18,],mB1r_h11_2[18,],mB1r_h12_2[18,],mB1r_h13_2[18,],
      mB1r_h20_2[18,],mB1r_h21_2[18,],mB1r_h22_2[18,],
      mBa1r_h10_1[18,],mBa1r_h11_1[18,],mBa1r_h12_1[18,],mBa1r_h13_1[18,],
      mBa1r_h20_1[18,],mBa1r_h21_1[18,],mBa1r_h22_1[18,],
      mBa1r_h10_2[18,],mBa1r_h11_2[18,],mBa1r_h12_2[18,],mBa1r_h13_2[18,],
      mBa1r_h20_2[18,],mBa1r_h21_2[18,],mBa1r_h22_2[18,],
      mBb1r_h10_1[18,],mBb1r_h11_1[18,],mBb1r_h12_1[18,],mBb1r_h13_1[18,],
      mBb1r_h20_1[18,],mBb1r_h21_1[18,],mBb1r_h22_1[18,],
      mBb1r_h10_2[18,],mBb1r_h11_2[18,],mBb1r_h12_2[18,],mBb1r_h13_2[18,],
      mBb1r_h20_2[18,],mBb1r_h21_2[18,],mBb1r_h22_2[18,],
      mBc1r_h10_1[18,],mBc1r_h11_1[18,],mBc1r_h12_1[18,],mBc1r_h13_1[18,],
      mBc1r_h20_1[18,],mBc1r_h21_1[18,],mBc1r_h22_1[18,],
      mBc1r_h10_2[18,],mBc1r_h11_2[18,],mBc1r_h12_2[18,],mBc1r_h13_2[18,],
      mBc1r_h20_2[18,],mBc1r_h21_2[18,],mBc1r_h22_2[18,]))
htest0$int <- factor(htest0$int, levels=unique(htest0$int))
htest0$dv <- factor(htest0$dv, levels=unique(htest0$dv))
htest0$cp <- factor(htest0$cp, levels=rev(unique(htest0$cp)))
htest0$lo95 <- htest0$Estimate - qnorm(0.975)*htest0$Std..Error
htest0$up95 <- htest0$Estimate + qnorm(0.975)*htest0$Std..Error
htest0$lo90 <- htest0$Estimate - qnorm(0.95)*htest0$Std..Error
htest0$up90 <- htest0$Estimate + qnorm(0.95)*htest0$Std..Error

#' 
#' ### 実験情報刺激効果と収入・イデオロギーの交差項係数による仮説検証（図A9）
#'   

p <- ggplot(subset(htest0, int%in%c("世帯収入","自己申告イデオロギー")), 
            aes(x=cp, y=Estimate)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lo95,ymax=up95, color=int, alpha=int), 
                width=0.25, position = position_dodge(width=-0.7)) + 
  geom_errorbar(aes(ymin=lo90,ymax=up90, color=int, alpha=int), 
                width=0, size=1.5, position = position_dodge(width=-0.7)) + 
  geom_point(aes(color=int, shape=int, alpha=int), size=3, 
             position = position_dodge(width=-0.7)) +
  facet_grid(h~dv, scales = "free_y", space = "free_y") + 
  coord_flip() + 
  scale_color_manual(name="収入・イデオロギー", values=rep("black",2)) + 
  # scale_color_brewer(name="収入・イデオロギー", type="qual", palette=2) + 
  scale_shape_discrete(name="収入・イデオロギー") + 
  scale_alpha_manual(name = "収入・イデオロギー", values=c(1,rep(0.3,1))) + 
  labs(x=NULL, y="実験群比較変数と収入・イデオロギーの交差項係数\n（従属変数は理想消費税率の平方根）", 
       caption="分析の詳細は回帰表を参照。太線は90%信頼区間、細線は95%信頼区間を示している。",
       subtitle = "対象商品") + 
  theme_bw() + theme(legend.position="bottom",
                     plot.subtitle = element_text(hjust=0.5),
                     strip.placement = "outside")

#+ fig.width=8, fig.height=6, dev="png", dpi=300, echo=FALSE
p

#+ eval=FALSE
ggsave("../out_v4/htest0_v1.png", p, width=8, height=6)

#' 
#' ### 実験情報刺激効果と収入・イデオロギーの交差項係数による仮説検証（その他のイデオロギー指標を含む）（図A12）
#'   

p <- ggplot(htest0, 
            aes(x=cp, y=Estimate)) + 
  geom_hline(aes(yintercept=0), linetype=2) + 
  geom_errorbar(aes(ymin=lo95,ymax=up95, color=int, alpha=int), 
                width=0.25, position = position_dodge(width=-0.7)) + 
  geom_errorbar(aes(ymin=lo90,ymax=up90, color=int, alpha=int), 
                width=0, size=1.5, position = position_dodge(width=-0.7)) + 
  geom_point(aes(color=int, shape=int, alpha=int), size=3, 
             position = position_dodge(width=-0.7)) +
  facet_grid(h~dv, scales = "free_y", space = "free_y") + 
  coord_flip() + 
  scale_color_manual(name="収入・イデオロギー", values=rep("black",5)) + 
  # scale_color_brewer(name="収入・イデオロギー", type="qual", palette=2) + 
  scale_shape_discrete(name="収入・イデオロギー") + 
  scale_alpha_manual(name = "収入・イデオロギー", values=c(1,rep(0.3,4))) + 
  labs(x=NULL, y="実験群比較変数と収入・イデオロギーの交差項係数\n（従属変数は理想消費税率の平方根）", 
       caption="分析の詳細は回帰表を参照。太線は90%信頼区間、細線は95%信頼区間を示している。",
       subtitle = "対象商品") + 
  theme_bw() + theme(legend.position="bottom",
                     plot.subtitle = element_text(hjust=0.5),
                     strip.placement = "outside") + 
  guides(color=guide_legend(nrow=2,byrow=TRUE),
         shape=guide_legend(nrow=2,byrow=TRUE),
         alpha=guide_legend(nrow=2,byrow=TRUE))

#+ warning=FALSE, fig.width=9, fig.height=7, dev="png", dpi=300, echo=FALSE
p

#+ warning=FALSE, eval=FALSE
ggsave("htest0_v2.png", p, width=9, height=7)

#'
#' ## 表のエクスポート
#'

#' 
#' ### 直接効果
#' 
screenreg(list(m0_1,m0_2), 
          override.se = list(coeftest(m0_1,vcovHC(m0_1,"HC2"))[,2],
                             coeftest(m0_2,vcovHC(m0_2,"HC2"))[,2]),
          override.pvalues = list(coeftest(m0_1,vcovHC(m0_1,"HC2"))[,4],
                                  coeftest(m0_2,vcovHC(m0_2,"HC2"))[,4]),
          symbol = "+",
          single.row=TRUE, digits = 3, stars = c(0.001,0.01,0.05,0.1),
          custom.coef.map = vnmap,
          custom.model.names = c("1:生活必需品","2:その他すべて"),
          caption = "理想消費税率の平方根に実験情報刺激が与える効果（重回帰分析）",
          caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
          label="basetab", dcolumn = TRUE, booktabs = TRUE, use.packages = FALSE,
          custom.note = "%stars. 最小二乗法による重回帰分析、ロバスト標準誤差使用．")
# texreg(list(m0_1,m0_2), 
#        override.se = list(coeftest(m0_1,vcovHC(m0_1,"HC2"))[,2],
#                           coeftest(m0_2,vcovHC(m0_2,"HC2"))[,2]),
#        override.pvalues = list(coeftest(m0_1,vcovHC(m0_1,"HC2"))[,4],
#                                coeftest(m0_2,vcovHC(m0_2,"HC2"))[,4]),
#        # file = "basetab.html", symbol = "&dagger;", 
#        file = "basetab.tex", symbol = "\\dagger", 
#        single.row=TRUE, digits = 3, stars = c(0.001,0.01,0.05,0.1),
#        custom.coef.map = vnmap,
#        custom.model.names = c("1:生活必需品","2:その他すべて"),
#        caption = "理想消費税率の平方根に実験情報刺激が与える効果（重回帰分析）",
#        caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
#        label="basetab", dcolumn = TRUE, booktabs = TRUE, use.packages = FALSE,
#        custom.note = "%stars. 最小二乗法による重回帰分析、ロバスト標準誤差使用．")
# tmp <- readLines("basetab.tex")
# tmp <- gsub("{dagger}","{\\dagger}", tmp, fixed=TRUE)
# writeLines(tmp, "basetab.tex", useBytes = TRUE)

#' 
#' ### 収入・自己申告イデオロギーによる条件付け（表A3）
#' 

## Full Table
screenreg(list(mA_1,mA_2,mB_1,mB_2), 
          override.se = list(coeftest(mA_1,vcovHC(mA_1,"HC2"))[,2],
                             coeftest(mA_2,vcovHC(mA_2,"HC2"))[,2],
                             coeftest(mB_1,vcovHC(mB_1,"HC2"))[,2],
                             coeftest(mB_2,vcovHC(mB_2,"HC2"))[,2]),
          override.pvalues = list(coeftest(mA_1,vcovHC(mA_1,"HC2"))[,4],
                                  coeftest(mA_2,vcovHC(mA_2,"HC2"))[,4],
                                  coeftest(mB_1,vcovHC(mB_1,"HC2"))[,4],
                                  coeftest(mB_2,vcovHC(mB_2,"HC2"))[,4]),
          symbol = "+",
          single.row=FALSE, digits = 3, stars = c(0.001,0.01,0.05,0.1),
          custom.coef.map = vnmap,
          custom.model.names = c("1:生活必需品","2:その他すべて",
                                 "3:生活必需品","4:その他すべて"),
          custom.header = list("世帯収入" = 1:2, "自己申告イデオロギー" = 3:4), 
          caption = "理想消費税率の平方根に実験情報刺激が与える効果と収入・イデオロギー（重回帰分析）",
          caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
          label="maintab_full", dcolumn = TRUE, booktabs = TRUE, use.packages = FALSE,
          custom.note = "%stars. 最小二乗法による重回帰分析、ロバスト標準誤差使用．")
# texreg(list(mA_1,mA_2,mB_1,mB_2), 
#        override.se = list(coeftest(mA_1,vcovHC(mA_1,"HC2"))[,2],
#                           coeftest(mA_2,vcovHC(mA_2,"HC2"))[,2],
#                           coeftest(mB_1,vcovHC(mB_1,"HC2"))[,2],
#                           coeftest(mB_2,vcovHC(mB_2,"HC2"))[,2]),
#        override.pvalues = list(coeftest(mA_1,vcovHC(mA_1,"HC2"))[,4],
#                                coeftest(mA_2,vcovHC(mA_2,"HC2"))[,4],
#                                coeftest(mB_1,vcovHC(mB_1,"HC2"))[,4],
#                                coeftest(mB_2,vcovHC(mB_2,"HC2"))[,4]),
#        # file = "maintab_full.html", symbol = "&dagger;", 
#        file = "maintab_full.tex", symbol = "\\dagger", 
#        single.row=FALSE, digits = 3, stars = c(0.001,0.01,0.05,0.1),
#        custom.coef.map = vnmap,
#        custom.model.names = c("1:生活必需品","2:その他すべて",
#                               "3:生活必需品","4:その他すべて"),
#        custom.header = list("世帯収入" = 1:2, "自己申告イデオロギー" = 3:4), 
#        caption = "理想消費税率の平方根に実験情報刺激が与える効果と収入・イデオロギー（重回帰分析）",
#        caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
#        label="maintab_full", dcolumn = TRUE, booktabs = TRUE, use.packages = FALSE,
#        custom.note = "%stars. 最小二乗法による重回帰分析、ロバスト標準誤差使用．")
# tmp <- readLines("maintab_full.tex")
# tmp <- gsub("{dagger}","{\\dagger}", tmp, fixed=TRUE)
# writeLines(tmp, "maintab_full.tex", useBytes = TRUE)

#' 
#' ### 他のイデオロギーによる条件付け
#' 

## 他のイデオロギー（生活必需品）
screenreg(list(mBc_1,mBa_1,mBb_1), 
          override.se = list(coeftest(mBc_1,vcovHC(mBc_1,"HC2"))[,2],
                             coeftest(mBa_1,vcovHC(mBa_1,"HC2"))[,2],
                             coeftest(mBb_1,vcovHC(mBb_1,"HC2"))[,2]),
          override.pvalues = list(coeftest(mBc_1,vcovHC(mBc_1,"HC2"))[,4],
                                  coeftest(mBa_1,vcovHC(mBa_1,"HC2"))[,4],
                                  coeftest(mBb_1,vcovHC(mBb_1,"HC2"))[,4]),
          symbol = "+",
          single.row=TRUE, digits = 3, stars = c(0.001,0.01,0.05,0.1),
          custom.coef.map = vnmap2,
          custom.model.names = c("政党支持","外交安全保障","権利機会平等"),
          caption = "生活必需品の理想消費税率の平方根に実験情報刺激が与える効果とイデオロギー（重回帰分析）",
          caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
          label="idetab1", dcolumn = TRUE, booktabs = TRUE, use.packages = FALSE,
          custom.note = "%stars. 最小二乗法による重回帰分析、ロバスト標準誤差使用．")
# texreg(list(mBc_1,mBa_1,mBb_1), 
#        override.se = list(coeftest(mBc_1,vcovHC(mBc_1,"HC2"))[,2],
#                           coeftest(mBa_1,vcovHC(mBa_1,"HC2"))[,2],
#                           coeftest(mBb_1,vcovHC(mBb_1,"HC2"))[,2]),
#        override.pvalues = list(coeftest(mBc_1,vcovHC(mBc_1,"HC2"))[,4],
#                                coeftest(mBa_1,vcovHC(mBa_1,"HC2"))[,4],
#                                coeftest(mBb_1,vcovHC(mBb_1,"HC2"))[,4]),
#        # file = "idetab1.html", symbol = "&dagger;", 
#        file = "idetab1.tex", symbol = "\\dagger", 
#        single.row=TRUE, digits = 3, stars = c(0.001,0.01,0.05,0.1),
#        custom.coef.map = vnmap2,
#        custom.model.names = c("政党支持","外交安全保障","権利機会平等"),
#        caption = "生活必需品の理想消費税率の平方根に実験情報刺激が与える効果とイデオロギー（重回帰分析）",
#        caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
#        label="idetab1", dcolumn = TRUE, booktabs = TRUE, use.packages = FALSE,
#        custom.note = "%stars. 最小二乗法による重回帰分析、ロバスト標準誤差使用．")
# tmp <- readLines("idetab1.tex")
# tmp <- gsub("{dagger}","{\\dagger}", tmp, fixed=TRUE)
# writeLines(tmp, "idetab1.tex", useBytes = TRUE)

## 他のイデオロギー（その他の商品）
screenreg(list(mBc_2,mBa_2,mBb_2), 
          override.se = list(coeftest(mBc_2,vcovHC(mBc_2,"HC2"))[,2],
                             coeftest(mBa_2,vcovHC(mBa_2,"HC2"))[,2],
                             coeftest(mBb_2,vcovHC(mBb_2,"HC2"))[,2]),
          override.pvalues = list(coeftest(mBc_2,vcovHC(mBc_2,"HC2"))[,4],
                                  coeftest(mBa_2,vcovHC(mBa_2,"HC2"))[,4],
                                  coeftest(mBb_2,vcovHC(mBb_2,"HC2"))[,4]),
          symbol = "+",
          single.row=TRUE, digits = 3, stars = c(0.001,0.01,0.05,0.1),
          custom.coef.map = vnmap2,
          custom.model.names = c("政党支持","外交安全保障","権利機会平等"),
          caption = "その他全ての商品の理想消費税率の平方根に実験情報刺激が与える効果とイデオロギー（重回帰分析）",
          caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
          label="idetab2", dcolumn = TRUE, booktabs = TRUE, use.packages = FALSE,
          custom.note = "%stars. 最小二乗法による重回帰分析、ロバスト標準誤差使用．")
# texreg(list(mBc_2,mBa_2,mBb_2), 
#        override.se = list(coeftest(mBc_2,vcovHC(mBc_2,"HC2"))[,2],
#                           coeftest(mBa_2,vcovHC(mBa_2,"HC2"))[,2],
#                           coeftest(mBb_2,vcovHC(mBb_2,"HC2"))[,2]),
#        override.pvalues = list(coeftest(mBc_2,vcovHC(mBc_2,"HC2"))[,4],
#                                coeftest(mBa_2,vcovHC(mBa_2,"HC2"))[,4],
#                                coeftest(mBb_2,vcovHC(mBb_2,"HC2"))[,4]),
#        # file = "idetab2.html", symbol = "&dagger;", 
#        file = "idetab2.tex", symbol = "\\dagger", 
#        single.row=TRUE, digits = 3, stars = c(0.001,0.01,0.05,0.1),
#        custom.coef.map = vnmap2,
#        custom.model.names = c("政党支持","外交安全保障","権利機会平等"),
#        caption = "その他全ての商品の理想消費税率の平方根に実験情報刺激が与える効果とイデオロギー（重回帰分析）",
#        caption.above = TRUE, fontsize = "scriptsize", float.pos = "ht!",
#        label="idetab2", dcolumn = TRUE, booktabs = TRUE, use.packages = FALSE,
#        custom.note = "%stars. 最小二乗法による重回帰分析、ロバスト標準誤差使用．")
# tmp <- readLines("idetab2.tex")
# tmp <- gsub("{dagger}","{\\dagger}", tmp, fixed=TRUE)
# writeLines(tmp, "idetab2.tex", useBytes = TRUE)

#+ warning=FALSE, eval=FALSE, echo=FALSE
# Exporting PDF & HTML Output
# In R Studio
# rmarkdown::render('main_analysis_tax_v4_replication.R',
#                   rmarkdown::pdf_document(latex_engine="xelatex"),
#                   clean=TRUE, encoding = "UTF-8")
# rmarkdown::render('main_analysis_tax_v4.R', rmarkdown::github_document(toc=TRUE),
#                   clean=FALSE, encoding = "UTF-8")
# tmp <- list.files(getwd())
# tmp <- tmp[grep("\\.spin\\.R$|\\.spin\\.Rmd$|\\.spin\\.Rmdtmp$|\\.utf8\\.md$|\\.knit\\.md$",tmp)]
# for (i in 1:length(tmp)) file.remove(paste0(getwd(),"/",tmp[i]))

